1 /*
2
3 PhyML: a program that computes maximum likelihood phylogenies from
4 DNA or AA homologous sequences.
5
6 Copyright (C) Stephane Guindon. Oct 2003 onward.
7
8 All parts of the source except where indicated are distributed under
9 the GNU public licence. See http://www.opensource.org for details.
10
11 */
12
13 /*
14 ** spr.c: Routines for performing SPR moves on the tree.
15 **
16 ** Wim Hordijk Last modified: 28 August 2006
17 ** Stephane Guindon 2007
18 */
19
20 #include "spr.h"
21
22 #ifdef BEAGLE
23 #include "beagle_utils.h"
24 #endif
25
26 /*********************************************************/
27 /* Sort list of SPR move by putting the shallowest moves first */
Sort_Spr_List_Depth(t_tree * tree)28 void Sort_Spr_List_Depth(t_tree *tree)
29 {
30 int i,j;
31 t_spr *buff;
32
33 for(i=0;i<tree->n_moves-1;i++)
34 {
35 for(j=i+1;j<tree->n_moves;j++)
36 {
37 if(tree->spr_list_one_edge[j]->depth_path < tree->spr_list_one_edge[i]->depth_path)
38 {
39 buff = tree->spr_list_one_edge[i];
40 tree->spr_list_one_edge[i] = tree->spr_list_one_edge[j];
41 tree->spr_list_one_edge[j] = buff;
42 }
43 }
44 }
45 }
46
47 /*********************************************************/
48 /*********************************************************/
49 /* Sort list of SPR move by putting the more likely moves first */
Sort_Spr_List_LnL(t_spr ** list,int list_size,t_tree * tree)50 void Sort_Spr_List_LnL(t_spr **list, int list_size, t_tree *tree)
51 {
52 int i,j;
53 t_spr *buff;
54
55 for(i=0;i<list_size-1;i++)
56 {
57 for(j=i+1;j<list_size;j++)
58 {
59 if(list[j]->lnL > list[i]->lnL)
60 {
61 buff = list[i];
62 list[i] = list[j];
63 list[j] = buff;
64 }
65 }
66 }
67 }
68
69
70 /*********************************************************/
71 /*********************************************************/
72 /* Sort list of SPR move by putting the more parsimonious moves first */
Sort_Spr_List_Pars(t_tree * tree)73 void Sort_Spr_List_Pars(t_tree *tree)
74 {
75 int i,j;
76 t_spr *buff;
77
78 for(i=0;i<tree->size_spr_list_one_edge-1;i++)
79 {
80 for(j=i+1;j<tree->size_spr_list_one_edge;j++)
81 {
82 if(tree->spr_list_one_edge[j]->pars < tree->spr_list_one_edge[i]->pars)
83 {
84 buff = tree->spr_list_one_edge[i];
85 tree->spr_list_one_edge[i] = tree->spr_list_one_edge[j];
86 tree->spr_list_one_edge[j] = buff;
87 }
88 }
89 }
90 }
91
92
93 /*********************************************************/
94 /*********************************************************/
95 /*********************************************************/
96
Randomize_Spr_List(t_tree * tree)97 void Randomize_Spr_List(t_tree *tree)
98 {
99 int i,j;
100 t_spr *buff;
101
102 for(i=0;i<tree->n_moves;i++)
103 {
104 j = Rand_Int(0,tree->n_moves-1);
105 buff = tree->spr_list_one_edge[i];
106 tree->spr_list_one_edge[i] = tree->spr_list_one_edge[j];
107 tree->spr_list_one_edge[j] = buff;
108 }
109 }
110
111 /*********************************************************/
112
Spr(phydbl init_lnL,phydbl prop_spr,t_tree * tree)113 int Spr(phydbl init_lnL, phydbl prop_spr, t_tree *tree)
114 {
115 int i,br;
116 int *br_idx;
117 t_edge *b;
118
119 tree->mod->s_opt->n_improvements = 0;
120 tree->mod->s_opt->max_spr_depth = 0;
121 tree->mod->s_opt->max_rank_triple_move = 0;
122 tree->mod->s_opt->max_delta_lnL_spr_current = 0.0;
123
124 Reset_Spr_List(tree->spr_list_all_edge,tree->size_spr_list_all_edge);
125
126 br_idx = Permutate(2*tree->n_otu-3);
127
128 Set_Both_Sides(YES,tree);
129 Lk(NULL,tree);
130 tree->best_lnL = tree->c_lnL;
131
132 /* PhyML_Printf("\n. init: %f",tree->c_lnL); */
133
134 for(i=0;i<MAX(1,(int)((2*tree->n_otu-3)*prop_spr));++i)
135 {
136 br = br_idx[i];
137
138 if(!(br%10)) if(tree->io->print_json_trace == YES) JSON_Tree_Io(tree,tree->io->fp_out_json_trace);
139
140 b = tree->a_edges[br];
141
142 if(b->l->v > tree->mod->s_opt->l_min_spr)
143 {
144 Spr_Subtree(b,b->left,tree);
145 Spr_Subtree(b,b->rght,tree);
146 }
147 }
148
149
150 Free(br_idx);
151
152 if(tree->mod->s_opt->n_improvements == 0 && tree->mod->s_opt->spr_lnL == YES)
153 {
154 Optimize_Br_Len_Serie(2,tree);
155 tree->best_lnL = tree->c_lnL;
156
157 Sort_Spr_List_LnL(tree->spr_list_all_edge,tree->size_spr_list_all_edge,tree);
158 for(i=0;i<MIN(20,2*tree->n_otu-3);++i)
159 {
160 Try_One_Spr_Move_Full(tree->spr_list_all_edge[i],NO,tree);
161 if(tree->spr_list_all_edge[i]->lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move) break;
162 }
163 Sort_Spr_List_LnL(tree->spr_list_all_edge,tree->size_spr_list_all_edge,tree);
164 if(tree->spr_list_all_edge[0]->lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move)
165 Try_One_Spr_Move_Full(tree->spr_list_all_edge[0],YES,tree);
166 }
167
168 return tree->mod->s_opt->n_improvements;
169 }
170
171 /*********************************************************/
172
Spr_Pre_Order(t_node * a,t_node * d,t_edge * b,t_tree * tree)173 void Spr_Pre_Order(t_node *a, t_node *d, t_edge *b, t_tree *tree)
174 {
175 if(d->tax) return;
176 else
177 {
178 unsigned int i;
179
180
181 /* printf("\n. a: %d d: %d score: %d d1: %d d2: %d ",a->num,d->num,tree->c_pars); */
182
183 /* for(i=0;i<3;++i) */
184 /* { */
185 /* if(d->v[i] != a) */
186 /* { */
187 /* Spr_Subtree(d->b[i],d->v[i],tree); */
188 /* } */
189 /* } */
190
191 Spr_Subtree(b,a,tree);
192
193 for(i=0;i<3;++i)
194 {
195 if(d->v[i] != a)
196 {
197 Spr_Pre_Order(d,d->v[i],d->b[i],tree);
198 }
199 }
200
201 }
202 }
203
204 /*********************************************************/
205
Spr_Subtree(t_edge * b,t_node * link,t_tree * tree)206 void Spr_Subtree(t_edge *b, t_node *link, t_tree *tree)
207 {
208 int i;
209 int n_moves_pars, n_moves, min_pars, best_move_idx;
210 t_spr *best_pars_move;
211 t_edge *init_target, *dummy, *residual;
212
213 if(link->v[0] == NULL || link->v[1] == NULL || link->v[2] == NULL) return;
214
215 Reset_Spr_List(tree->spr_list_one_edge,tree->size_spr_list_one_edge);
216
217 best_move_idx = -1;
218 tree->n_moves = 0;
219
220 MIXT_Set_Pars_Thresh(tree);
221
222 if((link != b->left) && (link != b->rght))
223 {
224 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
225 Exit("\n");
226 }
227 else
228 {
229 tree->mod->s_opt->worst_lnL_spr = BIG;
230
231 if(!link->tax) Test_All_Spr_Targets(b,link,tree);
232
233 if(tree->n_moves)
234 {
235 n_moves_pars = MIN(tree->mod->s_opt->min_n_triple_moves,tree->n_moves);
236 n_moves = MIN(tree->mod->s_opt->min_n_triple_moves,tree->n_moves);
237
238 if(tree->mod->s_opt->spr_lnL == NO) n_moves = n_moves_pars;
239 n_moves = MAX(1,n_moves);
240
241 if(tree->mod->s_opt->spr_pars == YES)
242 {
243 min_pars = 1E+8;
244 best_pars_move = NULL;
245
246 for(i=0;i<n_moves;++i)
247 {
248 if(tree->spr_list_one_edge[i]->pars < min_pars)
249 {
250 best_pars_move = tree->spr_list_one_edge[i];
251 min_pars = tree->spr_list_one_edge[i]->pars;
252 }
253 }
254
255 assert(best_pars_move);
256
257 if(best_pars_move->pars < tree->best_pars)
258 {
259 Prune_Subtree(best_pars_move->n_link,best_pars_move->n_opp_to_link,&init_target,&residual,tree);
260 Graft_Subtree(best_pars_move->b_target,best_pars_move->n_link,NULL,residual,NULL,tree);
261 if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
262 {
263 Prune_Subtree(best_pars_move->n_link,best_pars_move->n_opp_to_link,&dummy,&residual,tree);
264 Graft_Subtree(init_target,best_pars_move->n_link,NULL,residual,NULL,tree);
265 Set_Both_Sides(YES,tree);
266 Pars(NULL,tree);
267 }
268 else
269 {
270 if(best_pars_move->depth_path > tree->mod->s_opt->max_spr_depth) tree->mod->s_opt->max_spr_depth = best_pars_move->depth_path;
271 Set_Both_Sides(YES,tree);
272 Pars(NULL,tree);
273 tree->best_pars = tree->c_pars;
274 if(tree->best_pars != best_pars_move->pars)
275 {
276 PhyML_Fprintf(stderr,"\n== best_pars = %d move_pars = %d",tree->best_pars,best_pars_move->pars);
277 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
278 Exit("\n");
279 }
280 tree->mod->s_opt->n_improvements++;
281 }
282 Reset_Spr_List(tree->spr_list_one_edge,tree->size_spr_list_one_edge);
283 }
284 else
285 {
286 Set_Both_Sides(YES,tree);
287 Pars(NULL,tree);
288 }
289 }
290 else
291 {
292 int apply_move = NO;
293 phydbl accept_prob,u;
294
295 if(tree->mod->s_opt->spr_lnL == YES)
296 {
297 Sort_Spr_List_LnL(tree->spr_list_one_edge,tree->size_spr_list_one_edge,tree);
298
299 if(tree->spr_list_one_edge[0]->lnL > tree->best_lnL)
300 {
301 best_move_idx = 0;
302 }
303 else if(tree->mod->s_opt->eval_list_regraft == YES)
304 {
305 best_move_idx = Evaluate_List_Of_Regraft_Pos_Triple(tree->spr_list_one_edge,n_moves,tree);
306 }
307 else
308 {
309 best_move_idx = -1;
310 }
311 }
312 else
313 {
314 best_move_idx = Evaluate_List_Of_Regraft_Pos_Triple(tree->spr_list_one_edge,n_moves,tree);
315 }
316
317 if(best_move_idx > -1)
318 {
319 if(Are_Equal(tree->annealing_temp,0.0,1.E-3) == NO)
320 {
321 accept_prob = exp((tree->spr_list_one_edge[best_move_idx]->lnL - tree->best_lnL)/tree->annealing_temp);
322 u = Uni();
323 if(!(u > accept_prob)) apply_move = YES;
324 }
325 else
326 {
327 if(tree->spr_list_one_edge[best_move_idx]->lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move)
328 apply_move = YES;
329 }
330 }
331
332 if((best_move_idx > -1) && (apply_move == YES))
333 {
334 Try_One_Spr_Move_Triple(tree->spr_list_one_edge[best_move_idx],tree);
335 }
336 else
337 {
338 Pars(NULL,tree);
339 }
340 }
341 }
342 Reset_Spr_List(tree->spr_list_one_edge,tree->size_spr_list_one_edge);
343 }
344 }
345
346 /*********************************************************/
347
Test_All_Spr_Targets(t_edge * b_pulled,t_node * n_link,t_tree * tree)348 int Test_All_Spr_Targets(t_edge *b_pulled, t_node *n_link, t_tree *tree)
349 {
350 t_node *n_opp_to_link,*n_v1,*n_v2;
351 t_edge *b_target,*b_residual;
352 int i,dir1,dir2;
353 scalar_dbl *init_l_v1, *init_l_v2, *init_l_pulled;
354 scalar_dbl *init_v_v1, *init_v_v2, *init_v_pulled;
355 int best_found;
356 phydbl init_lnL;
357
358 if(tree->mixt_tree != NULL)
359 {
360 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
361 Exit("\n");
362 }
363
364 init_lnL = tree->c_lnL;
365 b_target = b_residual = NULL;
366 n_opp_to_link = (n_link == b_pulled->rght)?(b_pulled->left):(b_pulled->rght);
367
368 init_l_pulled = Duplicate_Scalar_Dbl(b_pulled->l);
369 init_v_pulled = Duplicate_Scalar_Dbl(b_pulled->l_var);
370
371 dir1 = dir2 = -1;
372 for(i=0;i<3;i++)
373 if(n_link->v[i] != n_opp_to_link)
374 {
375 if(dir1<0) dir1 = i;
376 else dir2 = i;
377 }
378
379 assert(dir1 > -1);
380 assert(dir2 > -1);
381
382 if(n_link->v[dir1]->num < n_link->v[dir2]->num)
383 {
384 n_v1 = n_link->v[dir1];
385 n_v2 = n_link->v[dir2];
386 init_l_v1 = Duplicate_Scalar_Dbl(n_link->b[dir1]->l);
387 init_l_v2 = Duplicate_Scalar_Dbl(n_link->b[dir2]->l);
388 init_v_v1 = Duplicate_Scalar_Dbl(n_link->b[dir1]->l_var);
389 init_v_v2 = Duplicate_Scalar_Dbl(n_link->b[dir2]->l_var);
390 }
391 else
392 {
393 n_v1 = n_link->v[dir2];
394 n_v2 = n_link->v[dir1];
395 init_l_v1 = Duplicate_Scalar_Dbl(n_link->b[dir2]->l);
396 init_l_v2 = Duplicate_Scalar_Dbl(n_link->b[dir1]->l);
397 init_v_v1 = Duplicate_Scalar_Dbl(n_link->b[dir2]->l_var);
398 init_v_v2 = Duplicate_Scalar_Dbl(n_link->b[dir1]->l_var);
399 }
400
401 if(!(n_v1->tax && n_v2->tax)) /*! Pruning is meaningless otherwise */
402 {
403 Prune_Subtree(n_link,n_opp_to_link,&b_target,&b_residual,tree);
404
405 if(tree->mod->s_opt->spr_lnL == YES) Update_PMat_At_Given_Edge(b_target,tree);
406
407 for(i=0;i<tree->size_spr_list_one_edge;++i) tree->spr_list_one_edge[i]->path_prev = NULL;
408
409 tree->edge_list = NULL;
410 tree->node_list = NULL;
411 best_found = NO;
412 tree->depth_curr_path = 0;
413 tree->curr_path[0] = b_target->left;
414 Test_One_Spr_Target_Recur(b_target->rght,
415 b_target->left,
416 b_pulled,n_link,b_residual,b_target,&best_found,NULL,tree);
417
418 if(best_found == NO || tree->perform_spr_right_away == NO)
419 {
420 tree->depth_curr_path = 0;
421 tree->curr_path[0] = b_target->rght;
422 Test_One_Spr_Target_Recur(b_target->left,
423 b_target->rght,
424 b_pulled,n_link,b_residual,b_target,&best_found,NULL,tree);
425 }
426
427 Graft_Subtree(b_target,n_link,NULL,b_residual,NULL,tree);
428
429 if((n_link->v[dir1] != n_v1) || (n_link->v[dir2] != n_v2)) PhyML_Printf("\n== Warning: -- SWITCH NEEDED -- ! \n");
430
431 Copy_Scalar_Dbl(init_l_v1,n_link->b[dir1]->l);
432 Copy_Scalar_Dbl(init_v_v1,n_link->b[dir1]->l_var);
433
434 Copy_Scalar_Dbl(init_l_v2,n_link->b[dir2]->l);
435 Copy_Scalar_Dbl(init_v_v2,n_link->b[dir2]->l_var);
436
437 Copy_Scalar_Dbl(init_l_pulled,b_pulled->l);
438 Copy_Scalar_Dbl(init_v_pulled,b_pulled->l_var);
439
440 if(tree->mod->s_opt->spr_pars == NO)
441 {
442 Update_PMat_At_Given_Edge(n_link->b[dir1],tree);
443 Update_PMat_At_Given_Edge(n_link->b[dir2],tree);
444 Update_PMat_At_Given_Edge(b_pulled,tree);
445 }
446
447 if(tree->mod->s_opt->spr_pars == NO)
448 {
449 Update_Partial_Lk(tree,b_pulled, n_link);
450 Update_Partial_Lk(tree,b_target, n_link);
451 Update_Partial_Lk(tree,b_residual,n_link);
452 }
453 else
454 {
455 Update_Partial_Pars(tree,b_pulled, n_link);
456 Update_Partial_Pars(tree,b_target, n_link);
457 Update_Partial_Pars(tree,b_residual,n_link);
458 }
459
460 t_ll *e_ll = tree->edge_list->head;
461 t_ll *n_ll = tree->node_list->head;
462 t_edge *e;
463 t_node *n;
464 do
465 {
466 assert(e_ll);
467 assert(n_ll);
468
469 e = (t_edge *)e_ll->v;
470 n = (t_node *)n_ll->v;
471
472 /* printf("\n. update on edge %d node %d",e->num,n->num); fflush(NULL); */
473 if(tree->mod->s_opt->spr_lnL)
474 Update_Partial_Lk(tree,e,n);
475 else
476 Update_Partial_Pars(tree,e,n);
477
478 e_ll = e_ll->next;
479 n_ll = n_ll->next;
480 }
481 while(e_ll != NULL);
482
483
484
485 Free_Linked_List(tree->edge_list);
486 Free_Linked_List(tree->node_list);
487 }
488
489 tree->c_lnL = init_lnL;
490
491 Free_Scalar_Dbl(init_l_v1);
492 Free_Scalar_Dbl(init_l_v2);
493 Free_Scalar_Dbl(init_l_pulled);
494
495 Free_Scalar_Dbl(init_v_v1);
496 Free_Scalar_Dbl(init_v_v2);
497 Free_Scalar_Dbl(init_v_pulled);
498
499 return 0;
500 }
501
502 /*********************************************************/
503
Test_One_Spr_Target_Recur(t_node * a,t_node * d,t_edge * pulled,t_node * link,t_edge * residual,t_edge * init_target,int * best_found,t_spr * prev_move,t_tree * tree)504 void Test_One_Spr_Target_Recur(t_node *a, t_node *d, t_edge *pulled, t_node *link, t_edge *residual, t_edge *init_target, int *best_found, t_spr *prev_move, t_tree *tree)
505 {
506 unsigned int i;
507 t_spr *move,*next_move;
508
509 move = next_move = NULL;
510
511 if(*best_found == YES && tree->perform_spr_right_away == YES) return;
512
513 if(d->tax) return;
514 else
515 {
516 for(i=0;i<3;++i)
517 {
518 if(d->v[i] != a)
519 {
520
521 if(tree->mod->s_opt->spr_pars == NO)
522 Update_Partial_Lk(tree,d->b[i],d);
523 else
524 Update_Partial_Pars(tree,d->b[i],d);
525
526 /* printf("\n push edge %d node %d",d->b[i]->num,d->num); fflush(NULL); */
527 Push_Bottom_Linked_List(d->b[i],&tree->edge_list,NO);
528 Push_Bottom_Linked_List(d,&tree->node_list,NO);
529
530 tree->depth_curr_path++;
531
532 tree->curr_path[tree->depth_curr_path] = d->v[i];
533
534 if((tree->depth_curr_path <= tree->mod->s_opt->max_depth_path) &&
535 (tree->depth_curr_path >= tree->mod->s_opt->min_depth_path))
536 {
537
538 move = Test_One_Spr_Target(d->b[i],pulled,link,residual,init_target,d,tree);
539
540 move->path_prev = prev_move;
541
542 if((tree->mod->s_opt->spr_pars == NO && move->lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move) ||
543 (tree->mod->s_opt->spr_pars == YES && move->pars < tree->best_pars))
544 {
545 *best_found = YES;
546 }
547 }
548
549 bool go_to_next =
550 (tree->depth_curr_path < tree->mod->s_opt->max_depth_path &&
551 ((tree->mod->s_opt->spr_pars == NO && move && move->lnL > tree->best_lnL - tree->mod->s_opt->max_delta_lnL_spr) ||
552 tree->mod->s_opt->spr_pars == YES));
553
554 /* bool go_to_next = tree->depth_curr_path < tree->mod->s_opt->max_depth_path; */
555
556 if(go_to_next == YES) Test_One_Spr_Target_Recur(d,d->v[i],pulled,link,residual,init_target,best_found,move,tree);
557
558 tree->depth_curr_path--;
559 }
560 }
561 }
562 }
563
564 /*********************************************************/
565
Test_One_Spr_Target(t_edge * b_target,t_edge * b_arrow,t_node * n_link,t_edge * b_residual,t_edge * init_target,t_node * polarity,t_tree * tree)566 t_spr *Test_One_Spr_Target(t_edge *b_target, t_edge *b_arrow, t_node *n_link, t_edge *b_residual, t_edge *init_target, t_node *polarity, t_tree *tree)
567 {
568 scalar_dbl *init_target_l, *init_arrow_l, *init_residual_l;
569 scalar_dbl *init_target_v, *init_arrow_v, *init_residual_v;
570 int i,dir_v0,dir_v1,dir_v2;
571 scalar_dbl *l0,*l1,*l2;
572 scalar_dbl *v0,*v1,*v2;
573 t_node *n1,*n2;
574 phydbl init_lnL;
575 int init_pars;
576 t_spr *move;
577 unsigned int rk;
578
579 if(tree->mixt_tree != NULL)
580 {
581 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
582 Exit("\n");
583 }
584
585 tree->n_moves++;
586
587 init_lnL = tree->c_lnL;
588 init_pars = tree->c_pars;
589
590 move = tree->spr_list_one_edge[tree->size_spr_list_one_edge];
591
592 if(move->init_target_l == NULL)
593 {
594 move->init_target_l = Duplicate_Scalar_Dbl(init_target->l);
595 move->init_target_v = Duplicate_Scalar_Dbl(init_target->l_var);
596 }
597 else
598 {
599 Copy_Scalar_Dbl(init_target->l, move->init_target_l);
600 Copy_Scalar_Dbl(init_target->l_var,move->init_target_v);
601 }
602
603
604 // Save edge lengths so that they can be recovered in the end
605 init_target_l = Duplicate_Scalar_Dbl(b_target->l);
606 init_target_v = Duplicate_Scalar_Dbl(b_target->l_var);
607
608 init_arrow_l = Duplicate_Scalar_Dbl(b_arrow->l);
609 init_arrow_v = Duplicate_Scalar_Dbl(b_arrow->l_var);
610
611 init_residual_l = Duplicate_Scalar_Dbl(b_residual->l);
612 init_residual_v = Duplicate_Scalar_Dbl(b_residual->l_var);
613
614 Graft_Subtree(b_target,n_link,NULL,b_residual,NULL,tree);
615
616 if(tree->mod->s_opt->spr_lnL == YES)
617 {
618 Update_PMat_At_Given_Edge(b_target,tree);
619 Update_PMat_At_Given_Edge(b_residual,tree);
620 Update_Partial_Lk(tree,b_arrow,n_link);
621 Lk(b_arrow,tree);
622 }
623 else
624 {
625 Update_Partial_Pars(tree,b_arrow,n_link);
626 Pars(b_arrow,tree);
627 }
628
629 n1 = (b_residual->left == n_link)?(b_residual->rght):(b_residual->left);
630 n2 = (b_target->left == n_link)?(b_target->rght):(b_target->left);
631 dir_v1 = dir_v2 = dir_v0 = -1;
632 for(i=0;i<3;++i)
633 {
634 if(n_link->v[i] == n1) dir_v1 = i;
635 else if(n_link->v[i] == n2) dir_v2 = i;
636 else dir_v0 = i;
637 }
638
639 l0 = Duplicate_Scalar_Dbl(n_link->b[dir_v0]->l);
640 v0 = Duplicate_Scalar_Dbl(n_link->b[dir_v0]->l_var);
641
642 if(n_link->v[dir_v1]->num > n_link->v[dir_v2]->num)
643 {
644 l1 = Duplicate_Scalar_Dbl(n_link->b[dir_v2]->l);
645 v1 = Duplicate_Scalar_Dbl(n_link->b[dir_v2]->l_var);
646 l2 = Duplicate_Scalar_Dbl(n_link->b[dir_v1]->l);
647 v2 = Duplicate_Scalar_Dbl(n_link->b[dir_v1]->l_var);
648 }
649 else
650 {
651 l1 = Duplicate_Scalar_Dbl(n_link->b[dir_v1]->l);
652 v1 = Duplicate_Scalar_Dbl(n_link->b[dir_v1]->l_var);
653 l2 = Duplicate_Scalar_Dbl(n_link->b[dir_v2]->l);
654 v2 = Duplicate_Scalar_Dbl(n_link->b[dir_v2]->l_var);
655 }
656
657 for(i=0;i<tree->depth_curr_path+1;++i) move->path[i] = tree->curr_path[i];
658
659 if(move->l0 != NULL)
660 {
661 Free_Scalar_Dbl(move->l0);
662 Free_Scalar_Dbl(move->v0);
663 }
664
665 if(move->l1 != NULL)
666 {
667 Free_Scalar_Dbl(move->l1);
668 Free_Scalar_Dbl(move->v1);
669 }
670
671 if(move->l2 != NULL)
672 {
673 Free_Scalar_Dbl(move->l2);
674 Free_Scalar_Dbl(move->v2);
675 }
676
677 move->l0 = l0;
678 move->v0 = v0;
679
680 move->l1 = l1;
681 move->v1 = v1;
682
683 move->l2 = l2;
684 move->v2 = v2;
685
686 move->depth_path = tree->depth_curr_path;
687 move->pars = tree->c_pars;
688 move->lnL = tree->c_lnL;
689 move->b_target = b_target;
690 move->n_link = n_link;
691 move->b_opp_to_link = b_arrow;
692 move->b_init_target = init_target;
693 move->dist = b_target->topo_dist_btw_edges;
694 move->n_opp_to_link = (n_link==b_arrow->left)?(b_arrow->rght):(b_arrow->left);
695
696 rk = Include_One_Spr_To_List_Of_Spr(tree->spr_list_one_edge,tree->size_spr_list_one_edge,move,tree);
697 Include_One_Spr_To_List_Of_Spr(tree->spr_list_all_edge,tree->size_spr_list_all_edge,move,tree);
698
699 Prune_Subtree(n_link,
700 (n_link==b_arrow->left)?(b_arrow->rght):(b_arrow->left),
701 &b_target,
702 &b_residual,
703 tree);
704
705
706 Copy_Scalar_Dbl(init_target_l,b_target->l);
707 Copy_Scalar_Dbl(init_target_v,b_target->l_var);
708
709 Copy_Scalar_Dbl(init_arrow_l,b_arrow->l);
710 Copy_Scalar_Dbl(init_arrow_v,b_arrow->l_var);
711
712 Copy_Scalar_Dbl(init_residual_l,b_residual->l);
713 Copy_Scalar_Dbl(init_residual_v,b_residual->l_var);
714
715 if(tree->mod->s_opt->spr_lnL == YES) Update_PMat_At_Given_Edge(b_target,tree);
716
717 tree->c_lnL = init_lnL;
718 tree->c_pars = init_pars;
719
720 Free_Scalar_Dbl(init_target_l);
721 Free_Scalar_Dbl(init_arrow_l);
722 Free_Scalar_Dbl(init_residual_l);
723 Free_Scalar_Dbl(init_target_v);
724 Free_Scalar_Dbl(init_arrow_v);
725 Free_Scalar_Dbl(init_residual_v);
726
727 return tree->spr_list_one_edge[rk];
728 }
729
730 /*********************************************************/
731
Global_Spr_Search(t_tree * tree)732 void Global_Spr_Search(t_tree *tree)
733 {
734 unsigned int i,iter;
735 phydbl best_lnL;
736 t_tree *best_tree;
737 time_t t_cur;
738 phydbl mean_delta_lnL_spr,max_delta_lnL_spr,tune_l_mult;
739
740 unsigned int no_improv = 0;
741 unsigned int last_best_found = 0;
742 unsigned int hit_zero_improv = 0;
743 unsigned int freq = 1;
744 const unsigned int round_freq = 10;
745
746 best_lnL = UNLIKELY;
747 tree->verbose = (tree->verbose == VL0) ? VL0 : VL1;
748 best_tree = Duplicate_Tree(tree);
749
750 if(tree->io->print_json_trace == YES) JSON_Tree_Io(tree,tree->io->fp_out_json_trace);
751
752 if(tree->verbose > VL0 && tree->io->quiet == NO) PhyML_Printf("\n\n. Score of initial tree: %.2f",tree->c_lnL);
753
754
755 tree->mod->s_opt->min_diff_lk_move = 1.E-1;
756 tree->mod->s_opt->min_diff_lk_local = 1.E-1;
757 Round_Optimize(tree,1000);
758 best_lnL = tree->c_lnL;
759 Copy_Tree(tree,best_tree);
760
761
762 if(tree->verbose > VL0 && tree->io->quiet == NO) PhyML_Printf("\n\n. Starting first round of SPRs...\n");
763
764 tree->mod->s_opt->max_depth_path = MAX(40,1+(int)(tree->n_otu/4));
765 tree->mod->s_opt->max_delta_lnL_spr = 1000;
766 tree->mod->s_opt->l_min_spr = 0.0;
767 tree->mod->s_opt->spr_lnL = YES;
768 tree->mod->s_opt->spr_pars = NO;
769 tree->mod->s_opt->min_diff_lk_move = 1.E-0;
770 tree->mod->s_opt->min_diff_lk_local = 1.E-0;
771 tree->perform_spr_right_away = YES;
772 tree->mod->s_opt->eval_list_regraft = NO;
773 tree->mod->s_opt->max_delta_lnL_spr_current = 0.0;
774 tree->mod->s_opt->min_n_triple_moves = 1;
775 mean_delta_lnL_spr = 0.0;
776 max_delta_lnL_spr = 0.0;
777 hit_zero_improv = 0;
778 tune_l_mult = 0.01;
779 best_lnL = tree->c_lnL;
780
781 iter = 0;
782 do
783 {
784
785 if(!(iter%freq))
786 {
787 for(int i=0;i<2*tree->n_otu-3;++i) MIXT_Multiply_Scalar_Dbl(tree->a_edges[i]->l,Rgamma((phydbl)(tune_l_mult*iter+1),(phydbl)(1./(tune_l_mult*iter+1))));
788 for(int i=0;i<2*tree->n_otu-3;++i) Set_Scalar_Dbl_Min_Thresh(tree->mod->l_min,tree->a_edges[i]->l);
789 for(int i=0;i<2*tree->n_otu-3;++i) Set_Scalar_Dbl_Max_Thresh(tree->mod->l_max,tree->a_edges[i]->l);
790 freq--;
791 if(freq < 1) freq=1;
792 }
793
794 if(iter> 10) tree->mod->s_opt->l_min_spr = 1.E-3;
795
796 Spr(tree->c_lnL,1.0,tree);
797 Optimize_Br_Len_Serie(2,tree);
798
799 if(!(iter%round_freq) && iter > 0) Round_Optimize(tree,1000);
800
801 if(tree->verbose > VL0 && tree->io->quiet == NO)
802 {
803 time(&t_cur);
804 PhyML_Printf("\n\t%8ds | %3d | lnL=%12.1f | depth=%5d/%5d | improvements=%4d | delta_lnL=%7.1f/%7.1f %c",
805 (int)(t_cur-tree->t_beg),
806 iter+1,
807 tree->c_lnL,
808 tree->mod->s_opt->max_spr_depth,
809 tree->mod->s_opt->max_depth_path,
810 tree->mod->s_opt->n_improvements,
811 tree->mod->s_opt->max_delta_lnL_spr_current,
812 tree->mod->s_opt->max_delta_lnL_spr,
813 (tree->numerical_warning == YES) ? '!' : ' ');
814 }
815
816 if(tree->mod->s_opt->n_improvements > (int)(tree->n_otu/5))
817 {
818 tune_l_mult *= 1.2;
819 }
820
821 tree->mod->s_opt->max_depth_path = MAX(5,MAX(tree->mod->s_opt->max_spr_depth+6,(int)(0.6*tree->mod->s_opt->max_depth_path)));
822
823 if((iter%4) > 0 || iter == 0)
824 {
825 mean_delta_lnL_spr += tree->mod->s_opt->max_delta_lnL_spr_current;
826 if(tree->mod->s_opt->max_delta_lnL_spr_current > max_delta_lnL_spr) max_delta_lnL_spr = tree->mod->s_opt->max_delta_lnL_spr_current;
827 }
828 else if(iter > 0)
829 {
830 mean_delta_lnL_spr /= 4.0;
831 /* tree->mod->s_opt->max_delta_lnL_spr = MAX(50.,MIN(2.0*mean_delta_lnL_spr,0.5*tree->mod->s_opt->max_delta_lnL_spr)); */
832 tree->mod->s_opt->max_delta_lnL_spr = MAX(50.,2.*max_delta_lnL_spr);
833 mean_delta_lnL_spr = tree->mod->s_opt->max_delta_lnL_spr_current;
834 max_delta_lnL_spr = 0.0;
835 }
836
837
838 if(tree->c_lnL > best_lnL)
839 {
840 no_improv = 0;
841 best_lnL = tree->c_lnL;
842 Copy_Tree(tree,best_tree);
843 if(tree->verbose > VL0 && tree->io->quiet == NO) PhyML_Printf(" +");
844 if(tree->io->print_json_trace == YES) JSON_Tree_Io(tree,tree->io->fp_out_json_trace);
845 }
846
847 if(tree->mod->s_opt->n_improvements == 0)
848 {
849 hit_zero_improv++;
850 }
851
852 no_improv++;
853 iter++;
854 }
855 while(tree->mod->s_opt->n_improvements > 15 && no_improv < 10);
856
857
858 if(tree->verbose > VL0 && tree->io->quiet == NO) PhyML_Printf("\n\n. Second round of optimization...\n");
859
860
861 tree->mod->s_opt->max_depth_path = MAX(10,(int)(1.5*tree->mod->s_opt->max_depth_path));
862 tree->mod->s_opt->l_min_spr = 1.E-3;
863 tree->mod->s_opt->spr_lnL = YES;
864 tree->mod->s_opt->spr_pars = NO;
865 tree->mod->s_opt->min_diff_lk_move = 1.E-2;
866 tree->mod->s_opt->min_diff_lk_local = 1.E-2;
867 tree->mod->s_opt->apply_spr_right_away = YES;
868 tree->mod->s_opt->apply_spr = YES;
869 tree->mod->s_opt->eval_list_regraft = NO;
870 tree->mod->s_opt->max_delta_lnL_spr_current = 0.0;
871 tree->mod->s_opt->min_n_triple_moves = 1;
872 tree->mod->s_opt->deepest_path = 0;
873 mean_delta_lnL_spr = 0.0;
874 max_delta_lnL_spr = 0.0;
875 hit_zero_improv = 0;
876 no_improv = 0;
877
878 do
879 {
880 if(!(iter%freq))
881 {
882 for(int i=0;i<2*tree->n_otu-3;++i) MIXT_Multiply_Scalar_Dbl(tree->a_edges[i]->l,Rgamma((phydbl)(tune_l_mult*iter+1),(phydbl)(1./(tune_l_mult*iter+1))));
883 for(int i=0;i<2*tree->n_otu-3;++i) Set_Scalar_Dbl_Min_Thresh(tree->mod->l_min,tree->a_edges[i]->l);
884 for(int i=0;i<2*tree->n_otu-3;++i) Set_Scalar_Dbl_Max_Thresh(tree->mod->l_max,tree->a_edges[i]->l);
885 freq--;
886 if(freq < 1) freq=1;
887 }
888
889 Spr(tree->c_lnL,1.0,tree);
890 Optimize_Br_Len_Serie(2,tree);
891
892 if(!(iter%round_freq) && iter > 0) Round_Optimize(tree,1000);
893
894 if(tree->verbose > VL0 && tree->io->quiet == NO)
895 {
896 time(&t_cur);
897 PhyML_Printf("\n\t%8ds | %3d | lnL=%12.1f | depth=%5d/%5d | improvements=%4d | delta_lnL=%7.1f/%7.1f %c",
898 (int)(t_cur-tree->t_beg),
899 iter+1,
900 tree->c_lnL,
901 tree->mod->s_opt->max_spr_depth,
902 tree->mod->s_opt->max_depth_path,
903 tree->mod->s_opt->n_improvements,
904 tree->mod->s_opt->max_delta_lnL_spr_current,
905 tree->mod->s_opt->max_delta_lnL_spr,
906 (tree->numerical_warning == YES) ? '!' : ' ');
907 }
908
909 tree->mod->s_opt->max_depth_path = MAX(5,MAX(tree->mod->s_opt->max_spr_depth+4,(int)(0.8*tree->mod->s_opt->max_depth_path)));
910 tree->mod->s_opt->max_depth_path = MIN(20,tree->mod->s_opt->max_depth_path);
911
912 if((iter%4) > 0 || iter == 0)
913 {
914 mean_delta_lnL_spr += tree->mod->s_opt->max_delta_lnL_spr_current;
915 if(tree->mod->s_opt->max_delta_lnL_spr_current > max_delta_lnL_spr) max_delta_lnL_spr = tree->mod->s_opt->max_delta_lnL_spr_current;
916 }
917 else if(iter > 0)
918 {
919 mean_delta_lnL_spr /= 4.0;
920 tree->mod->s_opt->max_delta_lnL_spr = MAX(20.,2.*max_delta_lnL_spr);
921 mean_delta_lnL_spr = tree->mod->s_opt->max_delta_lnL_spr_current;
922 max_delta_lnL_spr = 0.0;
923 }
924
925
926 if(tree->c_lnL > best_lnL)
927 {
928 no_improv = 0;
929 best_lnL = tree->c_lnL;
930 Copy_Tree(tree,best_tree);
931 if(tree->verbose > VL0 && tree->io->quiet == NO) PhyML_Printf(" +");
932 if(tree->io->print_json_trace == YES) JSON_Tree_Io(tree,tree->io->fp_out_json_trace);
933 }
934
935 if(tree->mod->s_opt->n_improvements == 0)
936 {
937 hit_zero_improv++;
938 }
939
940 no_improv++;
941 iter++;
942 }
943 while(tree->mod->s_opt->n_improvements > 5 && no_improv < 10);
944
945 if(tree->verbose > VL0 && tree->io->quiet == NO) PhyML_Printf("\n\n. Third round of optimization...\n");
946 last_best_found = 0;
947
948 tree->mod->s_opt->max_depth_path = MAX(10,tree->mod->s_opt->max_depth_path);
949 tree->mod->s_opt->max_delta_lnL_spr = MAX(100.,tree->mod->s_opt->max_delta_lnL_spr);
950 tree->mod->s_opt->spr_lnL = YES;
951 tree->mod->s_opt->spr_pars = NO;
952 tree->mod->s_opt->l_min_spr = 1.E-4;
953 tree->mod->s_opt->min_diff_lk_move = 1.E-2;
954 tree->mod->s_opt->min_diff_lk_local = 1.E-2;
955 tree->mod->s_opt->apply_spr_right_away = YES;
956 tree->mod->s_opt->apply_spr = YES;
957 tree->mod->s_opt->eval_list_regraft = YES;
958 tree->mod->s_opt->max_delta_lnL_spr_current = 0.0;
959 tree->mod->s_opt->min_n_triple_moves = 5;
960 tree->mod->s_opt->max_rank_triple_move = 0;
961 tree->mod->s_opt->max_no_better_tree_found = 10;
962
963 do
964 {
965 Spr(tree->c_lnL,1.0,tree);
966 Optimize_Br_Len_Serie(2,tree);
967
968 if(!(iter%round_freq)) Round_Optimize(tree,1000);
969
970 if(tree->verbose > VL0 && tree->io->quiet == NO)
971 {
972 time(&t_cur);
973 PhyML_Printf("\n\t%8ds | %3d | lnL=%12.1f | depth=%5d/%5d | improvements=%4d | delta_lnL=%7.1f/%7.1f | triple moves=%4d %c",
974 (int)(t_cur-tree->t_beg),
975 iter+1,
976 tree->c_lnL,
977 tree->mod->s_opt->max_spr_depth,
978 tree->mod->s_opt->max_depth_path,
979 tree->mod->s_opt->n_improvements,
980 tree->mod->s_opt->max_delta_lnL_spr_current,
981 tree->mod->s_opt->max_delta_lnL_spr,
982 tree->mod->s_opt->min_n_triple_moves,
983 (tree->numerical_warning == YES) ? '!' : ' ');
984 }
985
986 tree->mod->s_opt->max_depth_path = MIN(30,MAX(5,MAX(2*tree->mod->s_opt->max_spr_depth,(int)(0.8*tree->mod->s_opt->max_depth_path))));
987 tree->mod->s_opt->max_delta_lnL_spr = MAX(100.,0.8*tree->mod->s_opt->max_delta_lnL_spr_current);
988 tree->mod->s_opt->min_n_triple_moves = MAX(5,2*tree->mod->s_opt->max_rank_triple_move);
989
990 if(tree->c_lnL > best_lnL)
991 {
992 last_best_found = 0;;
993 best_lnL = tree->c_lnL;
994 Copy_Tree(tree,best_tree);
995 if(tree->verbose > VL0 && tree->io->quiet == NO) PhyML_Printf(" +");
996 if(tree->io->print_json_trace == YES) JSON_Tree_Io(tree,tree->io->fp_out_json_trace);
997 }
998
999 last_best_found++;
1000 iter++;
1001 }
1002 while(tree->mod->s_opt->n_improvements > 0 && last_best_found <= tree->mod->s_opt->max_no_better_tree_found);
1003
1004 Copy_Tree(best_tree,tree);
1005
1006 if(tree->verbose > VL0 && tree->io->quiet == NO) PhyML_Printf("\n\n. Final optimisation steps...\n");
1007
1008 i = tree->verbose;
1009 tree->verbose = VL0;
1010 tree->mod->s_opt->min_diff_lk_move = 1.E-03;
1011 tree->mod->s_opt->min_diff_lk_local = 1.E-03;
1012 do
1013 {
1014 tree->mod->s_opt->fast_nni = NO;
1015 Round_Optimize(tree,10);
1016 if(!Check_NNI_Five_Branches(tree)) break;
1017 }
1018 while(1);
1019 tree->verbose = i;
1020
1021 Free_Tree(best_tree);
1022
1023 return;
1024 }
1025
1026 //////////////////////////////////////////////////////////////
1027 //////////////////////////////////////////////////////////////
1028
Speed_Spr(t_tree * tree,phydbl prop_spr,int max_cycles,phydbl delta_lnL)1029 void Speed_Spr(t_tree *tree, phydbl prop_spr, int max_cycles, phydbl delta_lnL)
1030 {
1031 int step,old_pars;
1032 phydbl old_lnL;
1033
1034 if(tree->lock_topo == YES)
1035 {
1036 PhyML_Fprintf(stderr,"\n== The tree topology is locked.");
1037 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1038 Exit("\n");
1039 }
1040
1041
1042 Set_Both_Sides(NO,tree);
1043 Pars(NULL,tree);
1044 if(tree->mod->s_opt->spr_pars == NO) Lk(NULL,tree);
1045 Record_Br_Len(tree);
1046
1047 tree->mod->s_opt->deepest_path = 0;
1048 tree->best_pars = tree->c_pars;
1049 tree->best_lnL = tree->c_lnL;
1050 old_lnL = tree->c_lnL;
1051 old_pars = tree->c_pars;
1052 step = 0;
1053 do
1054 {
1055 ++step;
1056
1057 old_lnL = tree->c_lnL;
1058 old_pars = tree->c_pars;
1059
1060 Set_Both_Sides(YES,tree);
1061 Pars(NULL,tree);
1062 if(tree->mod->s_opt->spr_pars == NO) Lk(NULL,tree);
1063 Spr(UNLIKELY,prop_spr,tree);
1064
1065 // Set maximum depth for future spr rounds to deepest spr found so far
1066 tree->mod->s_opt->max_depth_path = tree->mod->s_opt->max_spr_depth;
1067
1068 if(tree->mod->s_opt->spr_pars == NO)
1069 {
1070 if(tree->mod->s_opt->n_improvements > 0)
1071 {
1072 /* Optimise branch lengths */
1073 Optimize_Br_Len_Serie(2,tree);
1074 /* Print log-likelihood and parsimony scores */
1075 if(tree->verbose > VL2 && tree->io->quiet == NO) Print_Lk(tree,"[Branch lengths ]");
1076 }
1077 }
1078 else
1079 {
1080 if(tree->verbose > VL2 && tree->io->quiet == NO) Print_Pars(tree);
1081 }
1082
1083 Pars(NULL,tree);
1084
1085 if(tree->io->print_trace)
1086 {
1087 char *s = Write_Tree(tree);
1088 PhyML_Fprintf(tree->io->fp_out_trace,"[%f]%s\n",tree->c_lnL,s); fflush(tree->io->fp_out_trace);
1089 if((tree->io->print_site_lnl) && (!tree->mod->s_opt->spr_pars))
1090 {
1091 Print_Site_Lk(tree,tree->io->fp_out_lk);
1092 fflush(tree->io->fp_out_lk);
1093 }
1094 Free(s);
1095 }
1096
1097 if(tree->io->print_json_trace == YES) JSON_Tree_Io(tree,tree->io->fp_out_json_trace);
1098
1099 /* Record the current best log-likelihood and parsimony */
1100 if(tree->c_lnL > tree->best_lnL) tree->best_lnL = tree->c_lnL;
1101 if(tree->c_pars < tree->best_pars) tree->best_pars = tree->c_pars;
1102
1103 if(tree->mod->s_opt->spr_pars == NO)
1104 {
1105 if(tree->c_lnL < old_lnL-tree->mod->s_opt->min_diff_lk_local)
1106 {
1107 PhyML_Fprintf(stderr,"\n== old_lnL = %f c_lnL = %f",old_lnL,tree->c_lnL);
1108 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1109 Exit("");
1110 }
1111 }
1112 else
1113 {
1114 if(tree->c_pars > old_pars)
1115 {
1116 PhyML_Fprintf(stderr,"\n== old_pars = %d c_pars = %d",old_pars,tree->c_pars);
1117 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1118 Exit("");
1119 }
1120 }
1121
1122 /* Record the current best branch lengths */
1123 Record_Br_Len(tree);
1124
1125 /* Exit if no improvements after complete optimization */
1126 if(step+1 > max_cycles) break;
1127 if((tree->mod->s_opt->spr_pars == NO) && (fabs(old_lnL-tree->c_lnL) < delta_lnL)) break;
1128 if((tree->mod->s_opt->spr_pars == YES) && (fabs((phydbl)(old_pars-tree->c_pars)) < 1)) break;
1129 if(!tree->mod->s_opt->n_improvements) break;
1130 }
1131 while(1);
1132 }
1133
1134 /*********************************************************/
1135
Evaluate_List_Of_Regraft_Pos_Triple(t_spr ** spr_list,int list_size,t_tree * tree)1136 int Evaluate_List_Of_Regraft_Pos_Triple(t_spr **spr_list, int list_size, t_tree *tree)
1137 {
1138 t_spr *move;
1139 t_edge *init_target, *b_residual;
1140 int i,j,best_move,better_found;
1141 int dir_v0, dir_v1, dir_v2;
1142 scalar_dbl *recorded_l,*recorded_v;
1143 phydbl best_lnL,init_lnL;
1144 int recorded;
1145
1146 if(tree->mixt_tree != NULL)
1147 {
1148 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1149 Exit("\n");
1150 }
1151
1152 best_lnL = UNLIKELY;
1153 init_target = b_residual = NULL;
1154 best_move = -1;
1155 init_lnL = tree->c_lnL;
1156 recorded_v = recorded_l = NULL;
1157 better_found = NO;
1158
1159 if(list_size == 0)
1160 {
1161 PhyML_Fprintf(stderr,"\n== List size is 0 !");
1162 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1163 Exit("\n");
1164 }
1165
1166
1167 recorded = NO;
1168 for(i=0;i<list_size;i++)
1169 {
1170 move = spr_list[i];
1171
1172 if(!move)
1173 {
1174 PhyML_Fprintf(stderr,"\n== move is NULL\n");
1175 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1176 Exit("\n");
1177 }
1178
1179 if(move->b_target)
1180 {
1181 /* Record t_edge lengths */
1182 Record_Br_Len(tree);
1183
1184 /* Prune subtree */
1185 Prune_Subtree(move->n_link,
1186 move->n_opp_to_link,
1187 &init_target,
1188 &b_residual,
1189 tree);
1190
1191 if(recorded == NO)
1192 {
1193 /*! Rough optimisation of the branch length at prune site
1194 * We only need to perform this optimisation for the first
1195 * element of spr_list because the pruned subtree is the
1196 * same across all the elements of spr_list. It would not
1197 * be true in the general case
1198 */
1199 recorded = YES;
1200
1201 Br_Len_Opt(&(init_target->l->v),init_target,tree);
1202
1203 /*! Record branch length at prune site */
1204 if(recorded_l == NULL)
1205 {
1206 recorded_l = Duplicate_Scalar_Dbl(init_target->l);
1207 recorded_v = Duplicate_Scalar_Dbl(init_target->l_var);
1208 }
1209 else
1210 {
1211 Copy_Scalar_Dbl(init_target->l,recorded_l);
1212 Copy_Scalar_Dbl(init_target->l_var,recorded_v);
1213 }
1214
1215 Copy_Scalar_Dbl(recorded_l,move->init_target_l);
1216 Copy_Scalar_Dbl(recorded_v,move->init_target_v);
1217 }
1218 else
1219 {
1220 Copy_Scalar_Dbl(recorded_l,move->b_init_target->l);
1221 Copy_Scalar_Dbl(recorded_v,move->b_init_target->l_var);
1222
1223 Copy_Scalar_Dbl(recorded_l,move->init_target_l);
1224 Copy_Scalar_Dbl(recorded_v,move->init_target_v);
1225 }
1226
1227 /* Update the change proba matrix at prune position */
1228 Update_PMat_At_Given_Edge(init_target,tree);
1229
1230 /* Update conditional likelihoods along the path from the prune to
1231 the regraft position */
1232 MIXT_Set_Alias_Subpatt(YES,tree);
1233 Update_Partial_Lk_Along_A_Path(move->path,move->depth_path+1,tree);
1234 MIXT_Set_Alias_Subpatt(NO,tree);
1235
1236
1237 /* Regraft subtree */
1238 Graft_Subtree(move->b_target,move->n_link,NULL,b_residual,NULL,tree);
1239
1240 MIXT_Set_Alias_Subpatt(YES,tree);
1241 move->lnL = Triple_Dist(move->n_link,tree);
1242 MIXT_Set_Alias_Subpatt(NO,tree);
1243
1244 /* printf("\n. %d/%d init_lnL: %12G move->lnL= %12G best_lnL=%12G absolute_best=%12G", */
1245 /* i, */
1246 /* list_size, */
1247 /* init_lnL, */
1248 /* move->lnL, */
1249 /* best_lnL, */
1250 /* tree->best_lnL); */
1251
1252 /* Record updated branch lengths for this move */
1253 dir_v1 = dir_v2 = dir_v0 = -1;
1254 for(j=0;j<3;j++)
1255 {
1256 if(move->n_link->v[j] == move->n_opp_to_link) dir_v0 = j;
1257 else if(dir_v1 < 0) dir_v1 = j;
1258 else dir_v2 = j;
1259 }
1260
1261 Copy_Scalar_Dbl(move->n_link->b[dir_v0]->l, move->l0);
1262 Copy_Scalar_Dbl(move->n_link->b[dir_v0]->l_var,move->v0);
1263
1264 if(move->n_link->v[dir_v1]->num > move->n_link->v[dir_v2]->num)
1265 {
1266 Copy_Scalar_Dbl(move->n_link->b[dir_v2]->l, move->l1);
1267 Copy_Scalar_Dbl(move->n_link->b[dir_v2]->l_var,move->v1);
1268
1269 Copy_Scalar_Dbl(move->n_link->b[dir_v1]->l, move->l2);
1270 Copy_Scalar_Dbl(move->n_link->b[dir_v1]->l_var,move->v2);
1271 }
1272 else
1273 {
1274 Copy_Scalar_Dbl(move->n_link->b[dir_v1]->l, move->l1);
1275 Copy_Scalar_Dbl(move->n_link->b[dir_v1]->l_var,move->v1);
1276
1277 Copy_Scalar_Dbl(move->n_link->b[dir_v2]->l, move->l2);
1278 Copy_Scalar_Dbl(move->n_link->b[dir_v2]->l_var,move->v2);
1279 }
1280
1281 if(move->lnL > best_lnL)
1282 {
1283 best_lnL = move->lnL;
1284 best_move = i;
1285 }
1286
1287 /* Regraft the subtree at its original position */
1288 Prune_Subtree(move->n_link,
1289 move->n_opp_to_link,
1290 &move->b_target,
1291 &b_residual,
1292 tree);
1293
1294 Graft_Subtree(init_target,
1295 move->n_link,
1296 NULL,
1297 b_residual,
1298 NULL,
1299 tree);
1300
1301 /* Restore branch lengths */
1302 Restore_Br_Len(tree);
1303
1304 /* Update relevant change proba matrices */
1305 Update_PMat_At_Given_Edge(move->b_target,tree);
1306
1307 tree->c_lnL = init_lnL;
1308 }
1309
1310 /* PhyML_Printf("\n. [ %4d/%4d ] %f %f %s", */
1311 /* i,list_size,tree->best_lnL,move->lnL, */
1312 /* (move->lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move) ? "**" : ""); */
1313
1314
1315 /* Bail out as soon as you've found a true improvement */
1316 if(move->lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move)
1317 {
1318 better_found = YES;
1319 if(i > tree->mod->s_opt->max_rank_triple_move) tree->mod->s_opt->max_rank_triple_move = i;
1320 break;
1321 }
1322 }
1323
1324 /* PhyML_Printf("\n. max_improv = %f",max_improv); */
1325
1326
1327 if(better_found == NO)
1328 {
1329 MIXT_Set_Alias_Subpatt(YES,tree);
1330 for(i=0;i<list_size;i++)
1331 {
1332 move = spr_list[i];
1333 if(move->b_target)
1334 {
1335 for(j=0;j<3;++j) Update_PMat_At_Given_Edge(move->n_link->b[j],tree);
1336 for(j=0;j<3;++j) Update_Partial_Lk(tree,move->n_link->b[j],move->n_link);
1337
1338 /* TO DO : we don't need to update all these partial likelihoods here.
1339 Would need to record only those that were along the paths examined
1340 above */
1341
1342 for(j=0;j<3;++j)
1343 if(move->n_link->v[j] != move->n_opp_to_link)
1344 Pre_Order_Lk(move->n_link,move->n_link->v[j],tree);
1345
1346 break;
1347 }
1348 }
1349 MIXT_Set_Alias_Subpatt(NO,tree);
1350 }
1351
1352 #ifdef DEBUG
1353 if(best_move < 0 && list_size > 0)
1354 {
1355 PhyML_Printf("\n\n== Best_move < 0 !");
1356 PhyML_Printf("\n== List size = %d",list_size);
1357 PhyML_Printf("\n== Best lnL = %f",best_lnL);
1358 for(i=0;i<list_size;i++)
1359 {
1360 move = spr_list[i];
1361 PhyML_Printf("\n== move %p %p lnL: %f",move,move->b_target,move->lnL);
1362 }
1363
1364 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1365 Exit("\n");
1366 }
1367 #endif
1368
1369 Free_Scalar_Dbl(recorded_l);
1370 Free_Scalar_Dbl(recorded_v);
1371
1372 return best_move;
1373 }
1374
1375 /*********************************************************/
1376
Try_One_Spr_Move_Triple(t_spr * move,t_tree * tree)1377 int Try_One_Spr_Move_Triple(t_spr *move, t_tree *tree)
1378 {
1379 t_edge *init_target, *b_residual;
1380 int j;
1381 int dir_v0, dir_v1, dir_v2;
1382 int accept;
1383
1384 assert(move);
1385 if(move->n_link == NULL) return -1;
1386
1387 if(tree->mixt_tree != NULL)
1388 {
1389 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1390 Exit("\n");
1391 }
1392
1393 Record_Br_Len(tree);
1394
1395 Prune_Subtree(move->n_link,
1396 move->n_opp_to_link,
1397 &init_target,
1398 &b_residual,
1399 tree);
1400
1401 Copy_Scalar_Dbl(move->init_target_l,init_target->l);
1402 Copy_Scalar_Dbl(move->init_target_v,init_target->l_var);
1403
1404 Graft_Subtree(move->b_target,move->n_link,NULL,b_residual,NULL,tree);
1405
1406 dir_v1 = dir_v2 = dir_v0 = -1;
1407 for(j=0;j<3;j++)
1408 {
1409 if(move->n_link->v[j] == move->n_opp_to_link) dir_v0 = j;
1410 else if(dir_v1 < 0) dir_v1 = j;
1411 else dir_v2 = j;
1412 }
1413
1414 Copy_Scalar_Dbl(move->l0,move->n_link->b[dir_v0]->l);
1415 Copy_Scalar_Dbl(move->v0,move->n_link->b[dir_v0]->l_var);
1416
1417 if(move->n_link->v[dir_v1]->num > move->n_link->v[dir_v2]->num)
1418 {
1419 Copy_Scalar_Dbl(move->l1,move->n_link->b[dir_v2]->l);
1420 Copy_Scalar_Dbl(move->v1,move->n_link->b[dir_v2]->l_var);
1421
1422 Copy_Scalar_Dbl(move->l2,move->n_link->b[dir_v1]->l);
1423 Copy_Scalar_Dbl(move->v2,move->n_link->b[dir_v1]->l_var);
1424 }
1425 else
1426 {
1427 Copy_Scalar_Dbl(move->l1,move->n_link->b[dir_v1]->l);
1428 Copy_Scalar_Dbl(move->v1,move->n_link->b[dir_v1]->l_var);
1429
1430 Copy_Scalar_Dbl(move->l2,move->n_link->b[dir_v2]->l);
1431 Copy_Scalar_Dbl(move->v2,move->n_link->b[dir_v2]->l_var);
1432 }
1433
1434 accept = YES;
1435 if(!Check_Topo_Constraints(tree,tree->io->cstr_tree)) accept = NO;
1436
1437
1438 if(accept == YES) /* Apply the move */
1439 {
1440 time(&(tree->t_current));
1441 Pars(NULL,tree);
1442
1443 Update_PMat_At_Given_Edge(init_target,tree);
1444 for(int i=0;i<3;++i) Update_PMat_At_Given_Edge(move->n_link->b[i],tree);
1445 Post_Order_Lk(move->n_opp_to_link,move->n_link,tree);
1446 Pre_Order_Lk(move->n_opp_to_link,move->n_link,tree);
1447 Pre_Order_Lk(move->n_link,move->n_opp_to_link,tree);
1448 Lk(move->b_opp_to_link,tree);
1449
1450 if(fabs(tree->c_lnL - move->lnL) > tree->mod->s_opt->min_diff_lk_move)
1451 {
1452 PhyML_Fprintf(stderr,"\n== c_lnL = %f move_lnL = %f", tree->c_lnL,move->lnL);
1453 PhyML_Fprintf(stderr,"\n== %d l0=%G l1=%G l2=%G v0=%G v1=%G v2=%G",move->n_link->num,move->l0->v,move->l1->v,move->l2->v,move->v0->v,move->v1->v,move->v2->v);
1454 PhyML_Fprintf(stderr,"\n== Gamma MGF? %d",tree->io->mod->gamma_mgf_bl);
1455 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d.\n",__FILE__,__LINE__);
1456 Check_Lk_At_Given_Edge(YES,tree);
1457 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
1458 }
1459
1460 if(tree->verbose > VL2 && tree->io->quiet == NO)
1461 {
1462 Print_Lk_And_Pars(tree);
1463 PhyML_Printf(" [depth=%5d]",move->depth_path); fflush(NULL);
1464 }
1465
1466
1467 tree->mod->s_opt->n_improvements++;
1468
1469 t_spr *dum_move = move;
1470 phydbl delta = 0.0;
1471 while(dum_move)
1472 {
1473 delta = move->lnL - dum_move->lnL;
1474 if(delta > tree->mod->s_opt->max_delta_lnL_spr_current)
1475 tree->mod->s_opt->max_delta_lnL_spr_current = delta;
1476 dum_move = dum_move->path_prev;
1477 }
1478
1479
1480
1481 if(tree->c_lnL > tree->best_lnL) tree->best_lnL = tree->c_lnL;
1482 Record_Br_Len(tree);
1483
1484 if(move->depth_path > tree->mod->s_opt->deepest_path)
1485 tree->mod->s_opt->deepest_path = move->depth_path;
1486
1487 if(move->depth_path > tree->mod->s_opt->max_spr_depth) tree->mod->s_opt->max_spr_depth = move->depth_path;
1488
1489 return 1;
1490 }
1491 else // Go back to original topology
1492 {
1493 Prune_Subtree(move->n_link,
1494 move->n_opp_to_link,
1495 &move->b_target,
1496 &b_residual,
1497 tree);
1498
1499 Graft_Subtree(init_target,
1500 move->n_link,
1501 NULL,
1502 b_residual,
1503 NULL,
1504 tree);
1505
1506 Restore_Br_Len(tree);
1507
1508 return 0;
1509 }
1510 }
1511
1512 /*********************************************************/
1513
Try_One_Spr_Move_Full(t_spr * move,short int apply_move,t_tree * tree)1514 int Try_One_Spr_Move_Full(t_spr *move, short int apply_move, t_tree *tree)
1515 {
1516 t_edge *init_target, *b_residual;
1517 phydbl init_lnL;
1518
1519 assert(move);
1520 if(move->n_link == NULL) return -1;
1521 init_lnL = tree->c_lnL;
1522
1523 Record_Br_Len(tree);
1524
1525 Prune_Subtree(move->n_link,
1526 move->n_opp_to_link,
1527 &init_target,
1528 &b_residual,
1529 tree);
1530 Graft_Subtree(move->b_target,move->n_link,NULL,b_residual,NULL,tree);
1531
1532 Optimize_Br_Len_Serie(2,tree);
1533
1534 move->lnL = tree->c_lnL;
1535
1536 if(tree->c_lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move && apply_move == YES)
1537 {
1538 tree->best_lnL = tree->c_lnL;
1539 tree->mod->s_opt->n_improvements++;
1540 return 1;
1541 }
1542 else
1543 {
1544 Prune_Subtree(move->n_link,
1545 move->n_opp_to_link,
1546 &move->b_target,
1547 &b_residual,
1548 tree);
1549
1550 Graft_Subtree(init_target,
1551 move->n_link,
1552 NULL,
1553 b_residual,
1554 NULL,
1555 tree);
1556
1557 Restore_Br_Len(tree);
1558 tree->c_lnL = init_lnL;
1559 return 0;
1560 }
1561
1562 return -1;
1563 }
1564
1565 /*********************************************************/
1566
Include_One_Spr_To_List_Of_Spr(t_spr ** list,int list_size,t_spr * move,t_tree * tree)1567 unsigned int Include_One_Spr_To_List_Of_Spr(t_spr **list, int list_size, t_spr *move, t_tree *tree)
1568 {
1569 unsigned int i, rk;
1570 t_spr *buff_spr,*orig_move, *orig_move_list, *move_list;
1571 t_tree *orig_tree;
1572
1573 if(tree->mixt_tree != NULL)
1574 {
1575 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1576 Exit("\n");
1577 }
1578
1579 rk = 0;
1580
1581 if(((tree->mod->s_opt->spr_lnL == YES) && (move->lnL > list[list_size-1]->lnL)) ||
1582 ((tree->mod->s_opt->spr_lnL == NO) && (move->pars <= list[list_size-1]->pars)))
1583 {
1584 move_list = list[list_size-1];
1585
1586 /* printf("\n. Include move with lnL: %G to list %p",move->lnL,(void *)list); */
1587
1588 move_list->depth_path = move->depth_path;
1589 move_list->pars = move->pars;
1590 move_list->lnL = move->lnL;
1591 move_list->dist = move->dist;
1592
1593 if(move_list->l0 == NULL)
1594 {
1595 move_list->l0 = Duplicate_Scalar_Dbl(move->l0);
1596 move_list->v0 = Duplicate_Scalar_Dbl(move->v0);
1597 }
1598 else
1599 {
1600 Copy_Scalar_Dbl(move->l0,move_list->l0);
1601 Copy_Scalar_Dbl(move->v0,move_list->v0);
1602 }
1603
1604 if(move_list->l1 == NULL)
1605 {
1606 move_list->l1 = Duplicate_Scalar_Dbl(move->l1);
1607 move_list->v1 = Duplicate_Scalar_Dbl(move->v1);
1608 }
1609 else
1610 {
1611 Copy_Scalar_Dbl(move->l1,move_list->l1);
1612 Copy_Scalar_Dbl(move->v1,move_list->v1);
1613 }
1614
1615 if(move_list->l2 == NULL)
1616 {
1617 move_list->l2 = Duplicate_Scalar_Dbl(move->l2);
1618 move_list->v2 = Duplicate_Scalar_Dbl(move->v2);
1619 }
1620 else
1621 {
1622 Copy_Scalar_Dbl(move->l2,move_list->l2);
1623 Copy_Scalar_Dbl(move->v2,move_list->v2);
1624 }
1625
1626
1627 if(move_list->init_target_l == NULL)
1628 {
1629 move_list->init_target_l = Duplicate_Scalar_Dbl(move->init_target_l);
1630 move_list->init_target_v = Duplicate_Scalar_Dbl(move->init_target_v);
1631 }
1632 else
1633 {
1634 Copy_Scalar_Dbl(move->init_target_l,move_list->init_target_l);
1635 Copy_Scalar_Dbl(move->init_target_v,move_list->init_target_v);
1636 }
1637
1638 orig_move = move;
1639 orig_move_list = move_list;
1640 orig_tree = tree;
1641 do
1642 {
1643 move_list->b_target = move->b_target;
1644 move_list->n_link = move->n_link;
1645 move_list->n_opp_to_link = move->n_opp_to_link;
1646 move_list->b_opp_to_link = move->b_opp_to_link;
1647 move_list->b_init_target = move->b_init_target;
1648
1649 move = move->next;
1650 move_list = move_list->next;
1651 tree = tree->next;
1652 }
1653 while(tree);
1654
1655 move = orig_move;
1656 move_list = orig_move_list;
1657 tree = orig_tree;
1658
1659 for(i=0;i<move_list->depth_path+1;++i) move_list->path[i] = move->path[i];
1660
1661 for(i=list_size-1;i>0;i--)
1662 {
1663 if(((tree->mod->s_opt->spr_lnL == YES) && (list[i]->lnL > list[i-1]->lnL)) ||
1664 ((tree->mod->s_opt->spr_lnL == NO) && (list[i]->pars <= list[i-1]->pars)))
1665 {
1666
1667 orig_tree = tree;
1668 do
1669 {
1670 buff_spr = list[i-1];
1671 list[i-1] = list[i];
1672 list[i] = buff_spr;
1673
1674 if(tree->next) tree = tree->next;
1675 else tree = tree->next;
1676 }
1677 while(tree);
1678 tree = orig_tree;
1679
1680 }
1681 else
1682 {
1683 rk = i;
1684 break;
1685 }
1686 }
1687 }
1688 return rk;
1689 }
1690
1691 /*********************************************************/
1692
Random_Spr(int n_moves,t_tree * tree)1693 void Random_Spr(int n_moves, t_tree *tree)
1694 {
1695 int i;
1696 int br_pulled, br_target;
1697 t_spr *spr_struct;
1698 t_edge *target, *residual;
1699
1700 spr_struct = Make_One_Spr(tree);
1701 Init_One_Spr(spr_struct);
1702 target = residual = NULL;
1703
1704 for(i=0;i<n_moves;i++)
1705 {
1706 /* br_pulled = (int)((phydbl)rand()/RAND_MAX * (2*tree->n_otu-3-1)); */
1707 br_pulled = Rand_Int(0,2*tree->n_otu-3-1);
1708
1709 do
1710 {
1711 /* br_target = (int)((phydbl)rand()/RAND_MAX * (2*tree->n_otu-3-1)); */
1712 br_target = Rand_Int(0,2*tree->n_otu-3-1);
1713 }
1714 while(br_target == br_pulled);
1715
1716 spr_struct->n_link = tree->a_edges[br_pulled]->left;
1717 spr_struct->n_opp_to_link = tree->a_edges[br_pulled]->rght;
1718 spr_struct->b_opp_to_link = tree->a_edges[br_pulled];
1719 spr_struct->b_target = tree->a_edges[br_target];
1720 spr_struct->b_init_target = NULL;
1721
1722 if(!Check_Spr_Move_Validity(spr_struct,tree))
1723 {
1724 spr_struct->n_link = tree->a_edges[br_pulled]->rght;
1725 spr_struct->n_opp_to_link = tree->a_edges[br_pulled]->left;
1726 }
1727
1728 #ifdef DEBUG
1729 if(!Check_Spr_Move_Validity(spr_struct,tree))
1730 {
1731 Warn_And_Exit("\n== Could not find a valid move...\n");
1732 }
1733 #endif
1734
1735 if(spr_struct->n_link == spr_struct->b_target->left ||
1736 spr_struct->n_link == spr_struct->b_target->rght)
1737 {
1738 n_moves++;
1739 continue;
1740 }
1741
1742 Prune_Subtree(spr_struct->n_link,
1743 spr_struct->n_opp_to_link,
1744 &target,
1745 &residual,
1746 tree);
1747
1748 Graft_Subtree(spr_struct->b_target,
1749 spr_struct->n_link,
1750 NULL,
1751 residual,NULL,tree);
1752 }
1753 Free(spr_struct);
1754 }
1755
1756 /*********************************************************/
1757
Reset_Spr_List(t_spr ** list,int size)1758 void Reset_Spr_List(t_spr **list, int size)
1759 {
1760 int i;
1761
1762 for(i=0;i<size;++i)
1763 {
1764 list[i]->depth_path = 0;
1765 list[i]->pars = MAX_PARS;
1766 list[i]->lnL = UNLIKELY;
1767 list[i]->n_link = NULL;
1768 list[i]->n_opp_to_link = NULL;
1769 list[i]->b_target = NULL;
1770 }
1771 }
1772
1773 /*********************************************************/
1774
Check_Spr_Move_Validity(t_spr * this_spr_move,t_tree * tree)1775 int Check_Spr_Move_Validity(t_spr *this_spr_move, t_tree *tree)
1776 {
1777 int match;
1778
1779 match = 0;
1780 Found_In_Subtree(this_spr_move->n_link,
1781 this_spr_move->n_opp_to_link,
1782 this_spr_move->b_target->left,
1783 &match,
1784 tree);
1785
1786 if(match) return 0;
1787 else return 1;
1788 }
1789
1790 /*********************************************************/
1791
Spr_Pars(int threshold,int n_round_max,t_tree * tree)1792 void Spr_Pars(int threshold, int n_round_max, t_tree *tree)
1793 {
1794 int curr_pars,round;
1795
1796 if(tree->verbose > VL2 && tree->io->quiet == NO) PhyML_Printf("\n. Minimizing parsimony...\n");
1797
1798 tree->best_pars = 1E+8;
1799 tree->best_lnL = UNLIKELY;
1800 tree->mod->s_opt->spr_lnL = NO;
1801 tree->mod->s_opt->spr_pars = YES;
1802 curr_pars = tree->c_pars;
1803 tree->mod->s_opt->max_depth_path = tree->n_otu;
1804 round = 0;
1805 do
1806 {
1807 curr_pars = tree->c_pars;
1808 Speed_Spr(tree,1.0,1,0.0);
1809 }
1810 while(tree->mod->s_opt->n_improvements && fabs((phydbl)(tree->c_pars - curr_pars)) > threshold && round++ < n_round_max);
1811 }
1812
1813 //////////////////////////////////////////////////////////////
1814 //////////////////////////////////////////////////////////////
1815
Spr_Shuffle(t_tree * mixt_tree)1816 void Spr_Shuffle(t_tree *mixt_tree)
1817 {
1818 phydbl lk_old;
1819 int *orig_catg,n,n_iter;
1820 t_tree *tree,**tree_list;
1821
1822 if(mixt_tree->verbose > VL0) PhyML_Printf("\n\n. Refining the tree...\n");
1823
1824 /*! Get the number of classes in each mixture */
1825 orig_catg = MIXT_Get_Number_Of_Classes_In_All_Mixtures(mixt_tree);
1826
1827
1828 /*! Set the number of rate classes to (at most) 2.
1829 ! Propagate this to every mixture tree in the analysis
1830 */
1831 tree = mixt_tree;
1832 n = 0;
1833 do
1834 {
1835 #ifdef BEAGLE
1836 tree->b_inst = create_beagle_instance(tree, tree->io->quiet, tree->io);
1837 //Instead of capping the rate categories at 2, just
1838 //give the other categories 0 weight
1839 if(orig_catg[n] > 2) //should we even bother?
1840 {
1841 double cat_wghts[orig_catg[n]];
1842 //Give the first two categories equal weights
1843 cat_wghts[0] = 0.5;
1844 cat_wghts[1] = 0.5;
1845 int i;
1846 for(i=2;i<orig_catg[n];++i){
1847 cat_wghts[i] = 0.0;
1848 }
1849 int ret = beagleSetCategoryWeights(tree->b_inst,0,cat_wghts);
1850 if(ret<0) {fprintf(stderr, "beagleSetCategoryWeights() on instance %i failed:%i\n\n",tree->b_inst,ret);Exit(""); }
1851 tree->mod->optimizing_topology = true;
1852 }
1853 #endif
1854 tree->mod->ras->n_catg = MIN(2,orig_catg[n]);
1855 if(tree->mod->ras->invar == YES) tree->mod->ras->n_catg--;
1856 tree = tree->next_mixt;
1857 n++;
1858 }
1859 while(tree);
1860
1861
1862 /*! Make sure the number of trees in each mixture is at most 2
1863 */
1864 tree_list = MIXT_Record_All_Mixtures(mixt_tree);
1865 MIXT_Break_All_Mixtures(orig_catg,mixt_tree);
1866
1867 Set_Both_Sides(YES,mixt_tree);
1868 Lk(NULL,mixt_tree);
1869
1870
1871 mixt_tree->best_pars = 1E+8;
1872 mixt_tree->mod->s_opt->spr_lnL = NO;
1873 mixt_tree->mod->s_opt->spr_pars = NO;
1874 mixt_tree->mod->s_opt->quickdirty = NO;
1875 mixt_tree->best_lnL = mixt_tree->c_lnL;
1876 mixt_tree->mod->s_opt->max_depth_path = mixt_tree->n_otu;
1877 mixt_tree->mod->s_opt->min_diff_lk_move = 0.1;
1878 mixt_tree->annealing_temp = 0.0;
1879
1880 n_iter = 0;
1881 do
1882 {
1883 Set_Both_Sides(YES,mixt_tree);
1884 Lk(NULL,mixt_tree);
1885 Pars(NULL,mixt_tree);
1886 Record_Br_Len(mixt_tree);
1887
1888 mixt_tree->best_pars = mixt_tree->c_pars;
1889 mixt_tree->best_lnL = mixt_tree->c_lnL;
1890
1891 lk_old = mixt_tree->c_lnL;
1892 Spr(UNLIKELY,1.0,mixt_tree);
1893
1894 mixt_tree->annealing_temp -= 2.;
1895
1896 if(mixt_tree->annealing_temp < 0.0) mixt_tree->annealing_temp = 0.0;
1897
1898 if(mixt_tree->mod->s_opt->n_improvements < 5 ||
1899 mixt_tree->mod->s_opt->max_spr_depth < 2 ||
1900 FABS(lk_old-mixt_tree->c_lnL) < 5. ||
1901 ++n_iter > 10) break;
1902
1903 }
1904 while(1);
1905
1906 mixt_tree->annealing_temp = 0.0;
1907
1908 if(mixt_tree->verbose > VL0 && mixt_tree->io->quiet == NO)
1909 {
1910 PhyML_Printf("\n\n. End of refining stage...\n");
1911 }
1912
1913 /*! Go back to the original data structure, with potentially more
1914 ! than 2 trees per mixture
1915 */
1916 MIXT_Reconnect_All_Mixtures(tree_list,mixt_tree);
1917 Free(tree_list);
1918
1919 /*! Set the number of rate classes to their original values
1920 */
1921 tree = mixt_tree;
1922 n = 0;
1923 do
1924 {
1925 tree->mod->ras->n_catg = orig_catg[n];
1926 #ifdef BEAGLE
1927 tree->mod->optimizing_topology = false;
1928 //Reset the rate categories to their original weights
1929 if(orig_catg[n] > 2){
1930 update_beagle_ras(tree->mod);
1931 }
1932 #endif
1933 if(tree->mod->ras->invar == YES) tree->mod->ras->n_catg--;
1934 tree = tree->next_mixt;
1935 n++;
1936 }
1937 while(tree);
1938
1939 Free(orig_catg);
1940
1941 /*! Only the first two trees for each mixture have been modified so
1942 ! far -> need to update the other trees by copying the modified trees
1943 ! onto them.
1944 */
1945 tree = mixt_tree;
1946 do
1947 {
1948 if(tree != mixt_tree) Copy_Tree(mixt_tree,tree);
1949 tree = tree->next;
1950 }
1951 while(tree);
1952
1953 }
1954
1955 //////////////////////////////////////////////////////////////
1956 //////////////////////////////////////////////////////////////
1957
Spr_Random_Explore(t_tree * tree,phydbl anneal_temp,phydbl prop_spr,int do_rnd,int max_cycles)1958 void Spr_Random_Explore(t_tree *tree, phydbl anneal_temp, phydbl prop_spr, int do_rnd, int max_cycles)
1959 {
1960 int step,i,n_targets,n_rand,no_improvement;
1961 t_tree *best_tree;
1962 scalar_dbl **best_bl;
1963 t_node *rnd_node;
1964 t_edge *b_target,*b_residual,**target_list,*rnd_edge;
1965 phydbl true_best_lnL;
1966
1967 if(tree->lock_topo == YES)
1968 {
1969 PhyML_Fprintf(stderr,"\n== The tree topology is locked.");
1970 PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
1971 Exit("\n");
1972 }
1973
1974 Set_Both_Sides(NO,tree);
1975 Pars(NULL,tree);
1976 Lk(NULL,tree);
1977
1978 tree->mod->s_opt->max_depth_path = (int)(tree->n_otu/3);
1979 tree->mod->s_opt->min_diff_lk_move = 0.1;
1980 tree->mod->s_opt->spr_lnL = NO;
1981 tree->mod->s_opt->spr_pars = NO;
1982 tree->mod->s_opt->deepest_path = 0;
1983 tree->best_pars = tree->c_pars;
1984 step = 0;
1985 true_best_lnL = tree->c_lnL;
1986 best_tree = Make_Tree_From_Scratch(tree->n_otu,tree->data);
1987 best_bl = Copy_Br_Len(tree);
1988 target_list = (t_edge **)mCalloc(2*tree->n_otu-3,sizeof(t_edge *));
1989 n_targets = 0;
1990 no_improvement = 0;
1991 tree->annealing_temp = anneal_temp;
1992 Copy_Tree(tree,best_tree);
1993
1994 do
1995 {
1996
1997 if(do_rnd == YES)
1998 {
1999 n_rand = 0;
2000 do
2001 {
2002 rnd_node = tree->a_nodes[Rand_Int(tree->n_otu,2*tree->n_otu-3)];
2003 assert(rnd_node != tree->n_root && rnd_node->tax == NO);
2004
2005 rnd_edge = rnd_node->b[Rand_Int(0,2)];
2006
2007 Prune_Subtree(rnd_node,
2008 rnd_node == rnd_edge->left ? rnd_edge->rght : rnd_edge->left,
2009 &b_target,
2010 &b_residual,
2011 tree);
2012
2013 n_targets = 0;
2014 for(i=0;i<3;i++)
2015 if(b_target->left->v[i] != b_target->rght)
2016 Get_List_Of_Adjacent_Targets(b_target->left,b_target->left->v[i],NULL,&target_list,&n_targets,0,5);
2017
2018 for(i=0;i<3;i++)
2019 if(b_target->rght->v[i] != b_target->left)
2020 Get_List_Of_Adjacent_Targets(b_target->rght,b_target->rght->v[i],NULL,&target_list,&n_targets,0,5);
2021
2022 if(n_targets > 0) b_target = target_list[Rand_Int(0,n_targets-1)];
2023
2024 assert(b_target != NULL);
2025
2026 Graft_Subtree(b_target,rnd_node,NULL,b_residual,NULL,tree);
2027
2028 n_rand++;
2029 }
2030 while(n_rand != 1);
2031 }
2032
2033 Set_Both_Sides(YES,tree);
2034 Lk(NULL,tree);
2035 Pars(NULL,tree);
2036
2037 Print_Lk_And_Pars(tree);
2038
2039 if(tree->annealing_temp < 0.0) tree->annealing_temp = 0.0;
2040 if(prop_spr > 1.0) prop_spr = 1.0;
2041
2042 tree->best_lnL = tree->c_lnL;
2043 tree->best_pars = tree->c_pars;
2044 Spr(UNLIKELY,prop_spr,tree);
2045
2046 tree->annealing_temp -= 0.5;
2047 prop_spr+=0.2;
2048
2049 Optimiz_All_Free_Param(tree,(tree->io->quiet == YES)?(0):(tree->verbose > VL0));
2050 Optimize_Br_Len_Serie(2,tree);
2051
2052 if(tree->io->print_trace)
2053 {
2054 char *s = Write_Tree(tree);
2055 PhyML_Fprintf(tree->io->fp_out_trace,"[%f]%s\n",tree->c_lnL,s); fflush(tree->io->fp_out_trace);
2056 if((tree->io->print_site_lnl) && (!tree->mod->s_opt->spr_pars))
2057 {
2058 Print_Site_Lk(tree,tree->io->fp_out_lk);
2059 fflush(tree->io->fp_out_lk);
2060 }
2061 Free(s);
2062 }
2063
2064 if(tree->io->print_json_trace == YES) JSON_Tree_Io(tree,tree->io->fp_out_json_trace);
2065
2066 /* Record the current best log-likelihood and parsimony */
2067 if(tree->c_lnL > true_best_lnL)
2068 {
2069 no_improvement = 0;
2070 true_best_lnL = tree->c_lnL;
2071 For(i,2*tree->n_otu-1) Free_Scalar_Dbl(best_bl[i]);
2072 Free(best_bl);
2073 best_bl = Copy_Br_Len(tree);
2074 Copy_Tree(tree,best_tree); /* Record tree topology, branch lengths and model parameters */
2075 }
2076 else
2077 {
2078 no_improvement++;
2079 }
2080
2081 Transfer_Br_Len_To_Tree(best_bl,tree);
2082 Copy_Tree(best_tree,tree);
2083
2084 }
2085 while(++step <= max_cycles && tree->mod->s_opt->n_improvements > 0 && tree->mod->s_opt->max_spr_depth > 1);
2086
2087 Free(target_list);
2088 }
2089
2090 //////////////////////////////////////////////////////////////
2091 //////////////////////////////////////////////////////////////
2092
2093
Prune_Regraft_Time_Tree(t_tree * tree)2094 void Prune_Regraft_Time_Tree(t_tree *tree)
2095 {
2096 phydbl u,ratio;
2097 phydbl t_min,t_max;
2098 phydbl cur_lnL_seq,new_lnL_seq;
2099 phydbl cur_lnL_time,new_lnL_time;
2100 phydbl new_t;
2101 int i,j,k,prune_idx,n_regraft_nd,regraft_idx,dir_prune;
2102 phydbl *times;
2103 int rnd_dir,dir_v1,dir_v2,keepon;
2104 t_node *prune,*prune_daughter,*new_regraft_nd,*cur_regraft_nd;
2105 t_ll *regraft_nd_list;
2106 t_edge *target, *ori_target, *residual,*regraft_edge;
2107 phydbl regraft_t_min,regraft_t_max;
2108
2109 times = tree->times->nd_t;
2110
2111 do
2112 {
2113 keepon = NO;
2114 for(i=tree->n_otu;i<2*tree->n_otu-2;++i) // for each internal node
2115 {
2116
2117 TIMES_Update_Node_Ordering(tree);
2118
2119 RATES_Record_Times(tree);
2120
2121 cur_lnL_seq = tree->c_lnL;
2122 new_lnL_seq = UNLIKELY;
2123 cur_lnL_time = tree->times->c_lnL_times;
2124 new_lnL_time = UNLIKELY;
2125
2126 regraft_edge = NULL;
2127 new_regraft_nd = NULL;
2128 cur_regraft_nd = NULL;
2129 new_t = 0.0;
2130
2131 // Prune node
2132 prune_idx = i;
2133 prune = tree->a_nodes[prune_idx];
2134 assert(prune && prune->tax == NO);
2135
2136
2137 // Select a daughter of prune node
2138 dir_v1 = dir_v2 = -1;
2139 for(j=0;j<3;++j)
2140 if(prune->v[j] != prune->anc && prune->b[j] != tree->e_root)
2141 {
2142 if(dir_v1 < 0) dir_v1 = j;
2143 else dir_v2 = j;
2144 }
2145
2146 u = Uni();
2147 if(u < 0.5) rnd_dir = dir_v1;
2148 else rnd_dir = dir_v2;
2149
2150 prune_daughter = prune->v[rnd_dir];
2151 cur_regraft_nd = prune->v[rnd_dir == dir_v1 ? dir_v2 : dir_v1];
2152
2153 if(prune == tree->n_root)
2154 {
2155 if(prune_daughter == prune->v[dir_v1] && prune->v[dir_v2]->tax == YES)
2156 {
2157 prune_daughter = prune->v[dir_v2];
2158 cur_regraft_nd = prune->v[dir_v1];
2159 }
2160
2161 if(prune_daughter == prune->v[dir_v2] && prune->v[dir_v1]->tax == YES)
2162 {
2163 prune_daughter = prune->v[dir_v1];
2164 cur_regraft_nd = prune->v[dir_v2];
2165 }
2166 }
2167
2168 assert(prune_daughter->anc == prune);
2169
2170 dir_prune = -1;
2171 for(j=0;j<3;j++)
2172 {
2173 if(prune_daughter->v[j] == prune || prune_daughter->b[j] == tree->e_root)
2174 {
2175 dir_prune = j;
2176 break;
2177 }
2178 }
2179 assert(dir_prune > -1);
2180
2181
2182 // Get the list of potential regraft nodes (oldest node on regraft edge)
2183 regraft_nd_list = DATE_List_Of_Regraft_Nodes(prune_daughter->v[dir_prune],prune_daughter,®raft_t_min,®raft_t_max,NO,tree);
2184 assert(regraft_nd_list);
2185 if(prune == tree->n_root) Push_Bottom_Linked_List(prune,®raft_nd_list,YES);
2186
2187
2188 // Number of regraft nodes
2189 n_regraft_nd = Linked_List_Len(regraft_nd_list);
2190
2191
2192 for(j=0;j<n_regraft_nd;j++)
2193 {
2194 // Randomly select one (uniform)
2195 regraft_idx = Rand_Int(0,n_regraft_nd-1);
2196 new_regraft_nd = Linked_List_Elem(regraft_idx,regraft_nd_list);
2197
2198 // Time of regraft node
2199 t_max = MIN(times[prune_daughter->num],times[new_regraft_nd->num]);
2200 if(new_regraft_nd == tree->n_root) t_min = 10.0*t_max;
2201 else t_min = times[new_regraft_nd->anc->num];
2202 t_min = MAX(t_min,regraft_t_min);
2203
2204 new_t = Uni()*(t_max-t_min) + t_min;
2205
2206
2207 // New age
2208 if(prune == tree->n_root || new_regraft_nd == tree->n_root)
2209 {
2210 if(prune == tree->n_root)
2211 {
2212 if(prune_daughter == tree->n_root->v[1]) times[tree->n_root->num] = times[tree->n_root->v[2]->num];
2213 else times[tree->n_root->num] = times[tree->n_root->v[1]->num];
2214 times[prune_daughter->v[dir_prune]->num] = new_t;
2215 }
2216 if(new_regraft_nd == tree->n_root)
2217 {
2218 times[prune_daughter->v[dir_prune]->num] = times[tree->n_root->num];
2219 times[tree->n_root->num] = new_t;
2220 }
2221 }
2222 else
2223 {
2224 times[prune->num] = new_t;
2225 }
2226
2227
2228 // Prune
2229 target = residual = NULL;
2230 Prune_Subtree(prune_daughter->v[dir_prune],
2231 prune_daughter,
2232 &target,&residual,tree);
2233 ori_target = target;
2234
2235
2236 // Regraft edge is the one sitting above regraft_nd
2237 if(new_regraft_nd == tree->n_root->v[1] ||
2238 new_regraft_nd == tree->n_root->v[2] ||
2239 new_regraft_nd == tree->n_root) regraft_edge = tree->e_root;
2240 else
2241 {
2242 for(k=0;k<3;k++) if(new_regraft_nd->v[k] == new_regraft_nd->anc) break;
2243 assert(k!=3);
2244 regraft_edge = new_regraft_nd->b[k];
2245 }
2246
2247 assert(regraft_edge);
2248
2249
2250 // Regraft
2251 Graft_Subtree(regraft_edge,
2252 prune_daughter->v[dir_prune],
2253 prune_daughter,
2254 residual,
2255 new_regraft_nd,tree);
2256
2257
2258 if(!TIMES_Check_Node_Height_Ordering(tree))
2259 {
2260 PhyML_Fprintf(stderr,"\n== prune[%d]->t:%.3f daughter[%d]->t:%.3f prune_anc[%d]->t:%.3f regraft[%d]->t:%.3f regraft_anc[%d]->t:%.3f [effective:%d] t_prior_min/max: [prune:[%.3f %.3f] regraft:[%.3f %.3f]] ",
2261 prune->num,
2262 times[prune->num],
2263 prune_daughter->num,
2264 times[prune_daughter->num],
2265 prune->anc ? prune->anc->num : -1,
2266 prune->anc ? times[prune->anc->num] : -1.,
2267 new_regraft_nd->num,
2268 times[new_regraft_nd->num],
2269 new_regraft_nd->anc ? new_regraft_nd->anc->num : -1,
2270 new_regraft_nd->anc ? times[new_regraft_nd->anc->num] : +1.,
2271 prune->num,
2272 tree->times->t_prior_min[prune->num],
2273 tree->times->t_prior_max[prune->num],
2274 tree->times->t_prior_min[new_regraft_nd->num],
2275 tree->times->t_prior_max[new_regraft_nd->num]);
2276 PhyML_Fprintf(stderr,"\n== root: %d %d %d",tree->n_root->num,tree->n_root->v[1]->num,tree->n_root->v[2]->num);
2277 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
2278 }
2279
2280 DATE_Assign_Primary_Calibration(tree);
2281 new_lnL_time = TIMES_Lk_Times(NO,tree);
2282
2283 if(new_lnL_time > UNLIKELY)
2284 {
2285 Set_Both_Sides(NO,tree);
2286 new_lnL_seq = Lk(NULL,tree);
2287 }
2288
2289 ratio = (new_lnL_seq - cur_lnL_seq);
2290
2291
2292 if(ratio < .0)
2293 {
2294 // Reject
2295 Prune_Subtree(prune_daughter->v[dir_prune],
2296 prune_daughter,
2297 &target,&residual,tree);
2298 Graft_Subtree(ori_target,
2299 prune_daughter->v[dir_prune],
2300 prune_daughter,residual,prune == tree->n_root ? tree->n_root : cur_regraft_nd,tree);
2301
2302 RATES_Reset_Times(tree);
2303 RATES_Update_Cur_Bl(tree);
2304 DATE_Assign_Primary_Calibration(tree);
2305 TIMES_Lk_Times(NO,tree);
2306
2307
2308 if(!(tree->times->c_lnL_times > UNLIKELY))
2309 {
2310 printf("\n== time prune: %f",times[prune->num]);
2311 printf("\n== time prune_daughter: %f",times[prune_daughter->num]);
2312 printf("\n== prune: %d prune_daughter: %d prune_daughter->v[dir_prune]: %d cur_regraft_nd: %d new_regraft_nd: %d",
2313 prune->num,
2314 prune_daughter->num,
2315 prune_daughter->v[dir_prune]->num,
2316 cur_regraft_nd->num,
2317 new_regraft_nd->num);
2318 TIMES_Lk_Times(YES,tree);
2319 fflush(NULL);
2320 }
2321 assert(tree->times->c_lnL_times > UNLIKELY);
2322
2323 tree->c_lnL = cur_lnL_seq;
2324 tree->times->c_lnL_times = cur_lnL_time;
2325 }
2326 else
2327 {
2328 PhyML_Printf("\n. Hill-climbing step :: subtree [%4d/%4d] target [%4d/%4d] lnl: %f delta: %f",
2329 i,2*tree->n_otu-2,
2330 j,n_regraft_nd,
2331 tree->c_lnL,
2332 ratio);
2333 if(ratio > 10.) keepon = YES;
2334 break;
2335 }
2336 }
2337 Free_Linked_List(regraft_nd_list);
2338 }
2339 }while(keepon == YES);
2340 }
2341
2342
2343 //////////////////////////////////////////////////////////////
2344 //////////////////////////////////////////////////////////////
2345 //////////////////////////////////////////////////////////////
2346 //////////////////////////////////////////////////////////////
2347 //////////////////////////////////////////////////////////////
2348 //////////////////////////////////////////////////////////////
2349 //////////////////////////////////////////////////////////////
2350 //////////////////////////////////////////////////////////////
2351 //////////////////////////////////////////////////////////////
2352 //////////////////////////////////////////////////////////////
2353
2354 // ** EOF: spr.c
2355