Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

spr.c

Visit:

Last change on this file was 19480, checked in by westram, 17 months ago
reintegrates 'gde' into 'trunk' adds: log:branches/gde@19460:19479
File size: 117.1 KB

Line
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
23	/*
24	** BIG: Some big number.
25	*/
26
27	/* #define BIG 1e05 */
28
29	/*
30	** Global vars.
31	**
32	** - cur_lk: The current likelihood of the tree.
33	** - subtree_dist: The average subtree distances matrix.
34	** - seq_dist: The sequence distance matrix.
35	** - optim_cand: Array for holding candidate moves for local and global branch
36	** length optimization.
37	** - rgrft_cand: Array for holding candidate regraft positions.
38	** - v_tmp: The central t_node of the temporary regraft structure for
39	** estimating changes in likelihood.
40	** - path: The path through the tree during the recursive tree length
41	** calculation.
42	** - sum_scale_tmp Array for temporarily storing scaling factors.
43	** - p_lk_tmp: Temporary partial likelihood storage.
44	** - e_brent: A temporary t_edge to use for estimating distances using Brent.
45
46	** - tree->mod->s_opt->wim_n_rgrft: Number of promising regraft positions to consider when
47	performing all improving SPR moves.
48	** - tree->mod->s_opt->wim_n_optim: Number of candidate moves on which to perform local branch
49	length optimization.
50	** - tree->mod->s_opt->wim_max_dist: Maximum regraft distance to consider.
51	** - tree->mod->s_opt->wim_n_globl: Number of candidates moves on which to perform global branch
52	length optimization.
53	** - tree->mod->s_opt->wim_n_best: Number of promising regraft positions to consider when
54	performing only the best SPR move.
55
56	** - nr_d_l: Total number of change in tree length calculations done.
57	** - nr_d_lk: Total number of change in likelihood calculations done.
58	** - nr_loc: Total number of local branch length optimizations done.
59	** - nr_glb: Total number of global branch length optimizations done.
60	*/
61
62	phydbl cur_lk, *subtree_dist, sum_scale_tmp, *p_lk_tmp;
63	matrix *seq_dist;
64	_move_ optim_cand, rgrft_cand;
65	t_node v_tmp=NULL, *path;
66	t_edge *e_brent=NULL;
67	int nr_d_L, nr_d_lk, nr_loc, nr_glb;
68
69	/*
70	** Init_SPR: Initialize the SPR algorithm: allocate memory and set variables.
71	**
72	** Parameters:
73	** - tree: The current tree to use for initialization.
74	*/
75
76	void Init_SPR (t_tree *tree)
77	{
78	int i, nr_nodes, nr_edges;
79	t_node u_0, u_1, *u_2;
80
81	/*
82	** Get the SPR parameter values.
83	*/
84	nr_edges = 2*tree->n_otu-3;
85
86	if(tree->mod->s_opt->wim_n_rgrft < 0) tree->mod->s_opt->wim_n_rgrft = 1 + nr_edges / 5;
87	if(tree->mod->s_opt->wim_n_globl < 0) tree->mod->s_opt->wim_n_globl = 1 + nr_edges / 10;
88	if(tree->mod->s_opt->wim_max_dist < 0) tree->mod->s_opt->wim_max_dist = 1 + nr_edges / 10;
89	if(tree->mod->s_opt->wim_n_optim < 0) tree->mod->s_opt->wim_n_optim = 100;
90	if(tree->mod->s_opt->wim_n_best < 0) tree->mod->s_opt->wim_n_best = tree->mod->s_opt->wim_n_rgrft; /* can't
91	* be
92	* anything else
93	*/
94
95
96	/*
97	** If it doesn't exist yet, create the temporary regraft structure:
98	** a central t_node with three edges and tip nodes adjacent to it.
99	*/
100	if (v_tmp == NULL)
101	{
102
103	v_tmp=Make_Node_Light(0);
104	v_tmp->tax = 0;
105	u_0=Make_Node_Light(1);
106	u_0->tax = 1;
107	u_1=Make_Node_Light(2);
108	u_1->tax = 1;
109	u_2=Make_Node_Light(3);
110	u_2->tax = 1;
111
112	v_tmp->v[0] = u_0;
113	v_tmp->v[1] = u_1;
114	v_tmp->v[2] = u_2;
115	u_0->v[0] = v_tmp;
116	u_1->v[0] = v_tmp;
117	u_2->v[0] = v_tmp;
118
119	t_edge *edge_0 = Make_Edge_Light (v_tmp, u_0, 0);
120	Make_Edge_Lk (edge_0, tree);
121	t_edge *edge_1 = Make_Edge_Light (v_tmp, u_1,1);
122	Make_Edge_Lk (edge_1, tree);
123	t_edge *edge_2 = Make_Edge_Light (v_tmp, u_2,2);
124	Make_Edge_Lk (edge_2, tree);
125
126
127	/* For(i,tree->data->crunch_len) */
128	/* { */
129	/* For(j,tree->mod->ras->n_catg) */
130	/* { */
131	/* Free(edge_0->p_lk_rght[i][j]); */
132	/* } */
133	/* Free(edge_0->p_lk_rght[i]); */
134	/* } */
135	Free(edge_0->p_lk_rght);
136
137	if(!edge_0->rght->tax) Free(edge_0->sum_scale_rght);
138
139
140	/* For(i,tree->data->crunch_len) */
141	/* { */
142	/* For(j,tree->mod->ras->n_catg) */
143	/* { */
144	/* Free(edge_1->p_lk_rght[i][j]); */
145	/* } */
146	/* Free(edge_1->p_lk_rght[i]); */
147	/* } */
148	Free(edge_1->p_lk_rght);
149
150	if(!edge_1->rght->tax) Free(edge_1->sum_scale_rght);
151
152
153	/* For(i,tree->data->crunch_len) */
154	/* { */
155	/* For(j,tree->mod->ras->n_catg) */
156	/* { */
157	/* Free(edge_2->p_lk_rght[i][j]); */
158	/* } */
159	/* Free(edge_2->p_lk_rght[i]); */
160	/* } */
161	Free(edge_2->p_lk_rght);
162
163	if(!edge_2->rght->tax) Free(edge_2->sum_scale_rght);
164	}
165
166	/*
167	** If it doesn't exist yet, create the temporary edge.
168	*/
169	if (e_brent == NULL)
170	{
171
172	u_1=Make_Node_Light(4);
173	u_1->tax = 1;
174	u_2=Make_Node_Light(5);
175	u_2->tax = 1;
176	u_1->v[0] = u_2;
177	u_2->v[0] = u_1;
178	t_edge *edge_4 = Make_Edge_Light (u_1, u_2, 3);
179	Make_Edge_Lk (edge_4, tree);
180	e_brent = u_1->b[0];
181
182
183	/* For(i,tree->data->crunch_len) */
184	/* { */
185	/* For(j,tree->mod->ras->n_catg) */
186	/* { */
187	/* Free(edge_4->p_lk_rght[i][j]); */
188	/* } */
189	/* Free(edge_4->p_lk_rght[i]); */
190	/* } */
191	Free(edge_4->p_lk_rght);
192
193	if(!edge_4->rght->tax) Free(edge_4->sum_scale_rght);
194
195
196	/* For(i,tree->data->crunch_len) */
197	/* { */
198	/* For(j,tree->mod->ras->n_catg) */
199	/* { */
200	/* Free(edge_4->p_lk_left[i][j]); */
201	/* } */
202	/* Free(edge_4->p_lk_left[i]); */
203	/* } */
204	Free(edge_4->p_lk_left);
205
206	if(!edge_4->left->tax) Free(edge_4->sum_scale_left);
207	}
208
209	/*
210	** Allocate memory for temporarily storing partial likelihoods and
211	** scaling factors.
212	*/
213	p_lk_tmp = (phydbl )mCalloc (tree->n_patterntree->mod->ras->n_catg*tree->mod->ns, sizeof (phydbl));
214	/* p_lk_tmp = (phydbl *)mCalloc (tree->n_pattern, sizeof (phydbl )); */
215	/* for (i = 0; i < tree->n_pattern; i++) */
216	/* { */
217	/* p_lk_tmp[i] = (phydbl *)mCalloc (tree->mod->ras->n_catg, sizeof (phydbl )); */
218	/* for (j = 0; j < tree->mod->ras->n_catg; j++) */
219	/* { */
220	/* p_lk_tmp[i][j] = (phydbl )mCalloc (tree->mod->ns, sizeof (phydbl)); /
221	/* } */
222	/* } */
223	sum_scale_tmp = (phydbl *)mCalloc (tree->n_pattern, sizeof (phydbl));
224
225	/*
226	** Allocate memory for storing the average subtree distances.
227	*/
228	nr_nodes = 2*tree->n_otu-2;
229	subtree_dist = (phydbl *)malloc (nr_nodes sizeof (phydbl *));
230	for (i = 0; i < nr_nodes; i++)
231	{
232	subtree_dist[i] = (phydbl )malloc (nr_nodes sizeof (phydbl));
233	}
234
235	/*
236	** Allocate memory for storing the candidate regraft positions and
237	** t_edge length optimization moves.
238	*/
239	rgrft_cand = (_move_ *)malloc (MAX(tree->mod->s_opt->wim_n_rgrft,tree->mod->s_opt->wim_n_best) sizeof (_move_ *));
240	for (i = 0; i < MAX(tree->mod->s_opt->wim_n_rgrft,tree->mod->s_opt->wim_n_best); i++)
241	{
242	rgrft_cand[i] = (_move_ *)malloc (sizeof (_move_));
243	rgrft_cand[i]->path = (t_node *)malloc ((tree->mod->s_opt->wim_max_dist+2) sizeof (t_node *));
244	}
245	optim_cand = (_move_ *)malloc (tree->mod->s_opt->wim_n_optim sizeof (_move_ *));
246	for (i = 0; i < tree->mod->s_opt->wim_n_optim; i++)
247	{
248	optim_cand[i] = (_move_ *)malloc (sizeof (_move_));
249	optim_cand[i]->path = (t_node *)malloc ((tree->mod->s_opt->wim_max_dist+2) sizeof (t_node *));
250	}
251	path = (t_node *)malloc ((tree->mod->s_opt->wim_max_dist+2) sizeof (t_node *));
252
253	if(!tree->mat)
254	{
255	seq_dist = ML_Dist (tree->data, tree->mod);
256	tree->mat = seq_dist;
257	}
258	else
259	seq_dist = tree->mat;
260
261	/*
262	** Set variables.
263	*/
264	nr_d_L = 0;
265	nr_d_lk = 0;
266	nr_loc = 0;
267	nr_glb = 0;
268	}
269
270
271	/*
272	** Clean_SPR: Free up the used memory.
273	**
274	** Parameters:
275	** - tree: The current tree.
276	*/
277
278	void Clean_SPR (t_tree *tree)
279	{
280	int i;
281
282	/*
283	** Clean up the temporary regraft structure.
284	*/
285	Free_Node (v_tmp->v[0]);
286	Free_Node (v_tmp->v[1]);
287	Free_Node (v_tmp->v[2]);
288	v_tmp->b[0]->p_lk_rght = NULL;
289	Free_Edge_Lk (v_tmp->b[0]);
290	Free_Edge (v_tmp->b[0]);
291	v_tmp->b[1]->p_lk_rght = NULL;
292	Free_Edge_Lk(v_tmp->b[1]);
293	Free_Edge(v_tmp->b[1]);
294	v_tmp->b[2]->p_lk_rght = NULL;
295	Free_Edge_Lk (v_tmp->b[2]);
296	Free_Edge (v_tmp->b[2]);
297	Free_Node (v_tmp);
298	v_tmp = NULL;
299
300	/*
301	** Clean up the temporary edge.
302	*/
303	Free_Node (e_brent->left);
304	Free_Node (e_brent->rght);
305	e_brent->p_lk_left = NULL;
306	e_brent->p_lk_rght = NULL;
307	Free_Edge_Lk (e_brent);
308	Free_Edge (e_brent);
309	e_brent = NULL;
310
311	/*
312	** Free the temporary partial likelihood and scaling memory.
313	*/
314	/* for (i = 0; i < tree->n_pattern; i++) */
315	/* { */
316	/* for (j = 0; j < tree->mod->ras->n_catg; j++) */
317	/* { */
318	/* free (p_lk_tmp[i][j]); */
319	/* } */
320	/* free (p_lk_tmp[i]); */
321	/* } */
322	free (p_lk_tmp);
323	free (sum_scale_tmp);
324
325	/*
326	** Free the subtree distance matrix.
327	*/
328	for (i = 0; i < 2*tree->n_otu - 2; i++)
329	{
330	free (subtree_dist[i]);
331	}
332	free (subtree_dist);
333
334	/*
335	** Free the arrays for storing the candidate regrafting positions and
336	** t_edge length optimization moves.
337	*/
338	for (i = 0; i < MAX(tree->mod->s_opt->wim_n_rgrft,tree->mod->s_opt->wim_n_best); i++)
339	{
340	free (rgrft_cand[i]->path);
341	free (rgrft_cand[i]);
342	}
343	free (rgrft_cand);
344	for (i = 0; i < tree->mod->s_opt->wim_n_optim; i++)
345	{
346	free (optim_cand[i]->path);
347	free (optim_cand[i]);
348	}
349	free (optim_cand);
350	free (path);
351
352
353	/*
354	** Print some statistics (for "research" purposes only).
355	*/
356	/* PhyML_Printf ("nr_d_L: %d\n", nr_d_L); */
357	/* PhyML_Printf ("nr_d_lk: %d\n", nr_d_lk); */
358	/* PhyML_Printf ("nr_loc: %d\n", nr_loc); */
359	/* PhyML_Printf ("nr_glb: %d\n", nr_glb); */
360	}
361
362
363	/*
364	** Optim_SPR: Optimize the tree using SPR moves.
365	**
366	** Parameters:
367	** - tree: The tree to optimize.
368	** - max_size: The maximum size (= number of taxa) of the subtrees to be
369	** pruned. If m=0 or m>ntax, all possible prunings will be
370	** considered.
371	** - method: The optimization method to use ("ALL" or "BEST").
372	*/
373
374	void Optim_SPR (t_tree *tree, int max_size, int method)
375	{
376	int nr_moves, improvement;
377
378	if(tree->mod->s_opt->print) PhyML_Printf("\n\n. Starting SPR moves...\n");
379
380	/*
381	** Calculate the current likelihood value.
382	*/
383	Set_Both_Sides(YES,tree);
384	cur_lk = Lk(NULL,tree);
385	time(&(tree->t_current));
386	if(tree->mod->s_opt->print) Print_Lk(tree,"topology");
387
388	/*
389	** Optimize all t_edge lengths and calculate the new likelihood value.
390	*/
391	/* PhyML_Printf("\n. Optimizing t_edge lengths."); */
392	Optimize_Br_Len_Serie (tree);
393	Set_Both_Sides(YES,tree);
394	cur_lk = Lk(NULL,tree);
395	time(&(tree->t_current));
396	if(tree->mod->s_opt->print) Print_Lk(tree,"topology");
397
398	/*
399	** While improvements were found, perform another round of SPR moves.
400	*/
401	nr_moves = 0;
402	improvement = 1;
403	while (improvement)
404	{
405	/*
406	** Perform one round of SPR moves.
407	*/
408	if (method == ALL)
409	{
410	improvement = Perform_SPR_Moves (tree, max_size);
411	}
412	else if (method == BEST)
413	{
414	improvement = Perform_Best_SPR (tree, max_size);
415	}
416	else if (method == ONE)
417	{
418	improvement = Perform_One_SPR (tree, max_size);
419	}
420	else
421	{
422	PhyML_Printf ("\n. Unknown SPR optimization method, bailing out...\n");
423	exit (1);
424	}
425
426	/* If an improvement was found, update statistics. */
427	if(improvement)
428	{
429	nr_moves++;
430	if((nr_moves == 1) \|\| (nr_moves % 4 == 1))
431	{
432	/*
433	** Optimize model parameters.
434	*/
435	Optimiz_All_Free_Param (tree,(tree->io->quiet)?(0):(tree->mod->s_opt->print));
436	Set_Both_Sides(YES,tree);
437	Lk(NULL,tree);
438	}
439	}
440
441	/* Beg SG 28 May 2007 */
442	if(method == BEST \|\| method == ONE) break;
443	/* Beg SG 28 May 2007 */
444	}
445
446	if(tree->mod->s_opt->print) PhyML_Printf ("\n\n. Number of SPR moves: %d\n", nr_moves);
447
448	/*
449	** Perform a last round of optimization steps (for t_edge lengths).
450	*/
451	Round_Optimize(tree,tree->data,ROUND_MAX);
452	Check_NNI_Five_Branches(tree);
453	}
454
455
456	/*
457	** Perform_SPR_Moves: Perform a round of SPR moves on the tree. Prune each subtree in
458	** turn and calculate the change in tree length for each candidate
459	** regraft position. Estimate change in likelihood for the most
460	** promising moves, and perform all moves that result in an
461	** improvement. If no improvements were found at all, try local edge
462	** length optimization. If still no improvement, try global edge
463	** length optimization.
464	**
465	** Parameters:
466	** - tree: The tree to perform the SPR moves on.
467	** - max_size: The maximum size (= number of taxa) of the subtrees to be
468	** pruned. If m=0 or m>ntax, all possible prunings will be
469	** considered.
470	**
471	** Returns:
472	** If the current tree could be improved: 1.
473	** Otherwise: 0.
474	*/
475
476	int Perform_SPR_Moves (t_tree *tree, int max_size)
477	{
478	int nr_edges, i, j, candidate, improvement;
479	t_node root, v_prune;
480	t_edge *e_prune;
481
482	/*
483	** Calculate the average subtree distances.
484	*/
485	root = tree->a_nodes[0];
486	PostOrder_v (tree, root->v[2], root->b[2]);
487
488	/*
489	** Initialize the array of optimization candidates.
490	*/
491	for (i = 0; i < tree->mod->s_opt->wim_n_optim; i++)
492	{
493	optim_cand[i]->delta_lk = -1.0*BIG;
494	optim_cand[i]->d_L = -1.0*BIG;
495	}
496
497	/*
498	** Try all possible SPR moves and perform the ones that give an improvement.
499	*/
500	nr_edges = 2*tree->n_otu - 3;
501	cur_lk = tree->c_lnL;
502	improvement = 0;
503
504
505	/* PhyML_Printf("\n >>>>>>>>>>>>>>>>>>"); */
506	/* PhyML_Printf("\n. cur_lk = %f %f",cur_lk,Lk(NULL,tree)); */
507
508	/* PhyML_Printf ("\n. Trying SPR moves"); */
509	/* PhyML_Printf ("\n. - calculating tree distances and estimating likelihoods"); */
510
511	for(i = 0; i < nr_edges; i++)
512	{
513	/*
514	** Get the next prune edge.
515	*/
516	e_prune = tree->a_edges[i];
517	/*
518	** Try right subtree if appropriate.
519	*/
520	if (!e_prune->left->tax)
521	{
522	/*
523	** Clear the regraft candidate list.
524	*/
525	for (j = 0; j < tree->mod->s_opt->wim_n_rgrft; j++)
526	{
527	rgrft_cand[j]->d_L = -1.0*BIG;
528	}
529	v_prune = e_prune->left;
530	/* if ((max_size == 0) \|\| (e_prune->num_tax_rght <= max_size)) */
531	{
532	/*
533	** Calculate changes in tree length, and estimate changes in likelihood for
534	** the most promising candidates. Perform moves that give an improvement.
535	*/
536	Calc_Tree_Length (e_prune, v_prune, tree);
537	if ((candidate = Est_Lk_Change (e_prune, v_prune, tree)) >= 0)
538	{
539	improvement = 1;
540	Make_Move (rgrft_cand[candidate],0,tree);
541	/* PhyML_Printf("\n. Make simple move"); */
542	/* PhyML_Printf("\n. lk after simple move = %f",Lk(tree)); */
543	}
544	}
545	}
546	/*
547	** Try left subtree if appropriate.
548	*/
549	if (!e_prune->rght->tax)
550	{
551	/*
552	** Clear the regraft candidate list.
553	*/
554	for (j = 0; j < tree->mod->s_opt->wim_n_rgrft; j++)
555	{
556	rgrft_cand[j]->d_L = -1.0*BIG;
557	}
558	v_prune = e_prune->rght;
559	/* if ((max_size == 0) \|\| (e_prune->num_tax_left <= max_size)) */
560	{
561	/*
562	** Calculate changes in tree length, and estimate changes in likelihood for
563	** the most promising candidates. Perform moves that give an improvement.
564	*/
565	Calc_Tree_Length (e_prune, v_prune, tree);
566	if ((candidate = Est_Lk_Change (e_prune, v_prune, tree)) >= 0)
567	{
568	improvement = 1;
569	Make_Move (rgrft_cand[candidate],0,tree);
570	/* PhyML_Printf("\n. Make simple move"); */
571	/* PhyML_Printf("\n. lk after simple move = %f",Lk(tree)); */
572	}
573	}
574	}
575	}
576
577	/*
578	** If there was no improvement at all, try local t_edge length optimization at the
579	** regraft position.
580	*/
581
582	/* PhyML_Printf("\n. before local = %f %f",tree->c_lnL,Lk(tree)); */
583
584	if (!improvement)
585	{
586	/* PhyML_Printf ("\n. - performing local t_edge length optimizations"); */
587	if ((candidate = Find_Optim_Local (tree)) >= 0)
588	{
589	/* PhyML_Printf("\n. make local move"); */
590	improvement = 1;
591	Make_Move (optim_cand[candidate],1,tree);
592	/* PhyML_Printf("\n. lk after local move = %f",Lk(tree)); */
593	}
594	}
595
596
597	/*
598	** If there was still no improvement, try global t_edge length optimization.
599	*/
600	/* PhyML_Printf("\n. before global = %f %f",tree->c_lnL,Lk(tree)); */
601
602	if (!improvement)
603	{
604	/* PhyML_Printf ("\n. - performing global t_edge length optimization"); */
605	if ((candidate = Find_Optim_Globl (tree)) >= 0)
606	{
607	/* PhyML_Printf("\n. make global move"); */
608	improvement = 1;
609	Make_Move (optim_cand[candidate],2,tree);
610	/* PhyML_Printf("\n. lk after global move = %f",Lk(tree)); */
611	}
612	}
613
614	/* PhyML_Printf("\n. after all = %f %f",tree->c_lnL,Lk(tree)); */
615
616	/*
617	** Optimize all t_edge lengths again to make sure we got an updated
618	** likelihood value.
619	*/
620	Set_Both_Sides(YES,tree);
621	cur_lk = Lk(NULL,tree);
622	root = tree->a_nodes[0];
623	Optimize_Br_Len_Serie (tree);
624	Set_Both_Sides(YES,tree);
625	cur_lk = Lk(NULL,tree);
626	time(&(tree->t_current));
627	if(tree->mod->s_opt->print) Print_Lk(tree,"topoLOGy");
628
629	/*
630	** Return the result.
631	*/
632	return (improvement);
633	}
634
635
636	/*
637	** Perform_Best_SPR: Perform the best SPR move on the tree. Prune each subtree in
638	** turn and calculate the change in tree length for each candidate
639	** regraft position. Estimate change in likelihood for the most
640	** promising regraft positions, and store the best one. Then choose
641	** the best candidate over all moves. If no improving move can be
642	** found, try local t_edge length optimization, and if necessary
643	** global t_edge length optimization.
644	**
645	** Parameters:
646	** - tree: The tree to perform the SPR moves on.
647	** - max_size: The maximum size (= number of taxa) of the subtrees to be
648	** pruned. If m=0 or m>ntax, all possible prunings will be
649	** considered.
650	**
651	** Returns:
652	** If an improving move could be performed: 1.
653	** Otherwise: 0.
654	*/
655
656	int Perform_Best_SPR (t_tree *tree, int max_size)
657	{
658	int nr_edges, i, j, candidate, improvement;
659	t_node root, v_prune;
660	t_edge *e_prune;
661
662	/*
663	** Calculate the average subtree distances.
664	*/
665	root = tree->a_nodes[0];
666	PostOrder_v (tree, root->v[2], root->b[2]);
667
668	/*
669	** Initialize the array of optimization candidates.
670	*/
671	for (i = 0; i < tree->mod->s_opt->wim_n_optim; i++)
672	{
673	optim_cand[i]->delta_lk = -1.0*BIG;
674	optim_cand[i]->d_L = -1.0*BIG;
675	}
676
677	/*
678	** Try all possible SPR moves and perform the best one.
679	*/
680	nr_edges = 2*tree->n_otu - 3;
681	cur_lk = tree->c_lnL;
682	improvement = 0;
683	/* PhyML_Printf ("\n. Trying SPR moves"); */
684	/* PhyML_Printf ("\n. -calculating tree distances and estimating likelihoods"); */
685	for (i = 0; i < nr_edges; i++)
686	{
687	/*
688	** Get the next prune edge.
689	*/
690	e_prune = tree->a_edges[i];
691	/*
692	** Try right subtree if appropriate.
693	*/
694	if (!e_prune->left->tax)
695	{
696	/*
697	** Clear the regraft candidate list.
698	*/
699	for (j = 0; j < tree->mod->s_opt->wim_n_best; j++)
700	{
701	rgrft_cand[j]->d_L = -1.0*BIG;
702	}
703	v_prune = e_prune->left;
704	/* if ((max_size == 0) \|\| (e_prune->num_tax_rght <= max_size)) */
705	{
706	/*
707	** Calculate changes in tree length, and estimate changes in likelihood for
708	** the most promising candidates. Store the best one in the optimization list.
709	*/
710	Calc_Tree_Length (e_prune, v_prune, tree);
711	candidate = Best_Lk_Change (e_prune, v_prune, tree);
712	}
713	}
714	/*
715	** Try left subtree if appropriate.
716	*/
717	if (!e_prune->rght->tax)
718	{
719	/*
720	** Clear the regraft candidate list.
721	*/
722	for (j = 0; j < tree->mod->s_opt->wim_n_rgrft; j++)
723	{
724	rgrft_cand[j]->d_L = -1.0*BIG;
725	}
726	v_prune = e_prune->rght;
727	/* if ((max_size == 0) \|\| (e_prune->num_tax_left <= max_size)) */
728	{
729	/*
730	** Calculate changes in tree length, and estimate changes in likelihood for
731	** the most promising candidates. Perform moves that give an improvement.
732	*/
733	Calc_Tree_Length (e_prune, v_prune, tree);
734	candidate = Best_Lk_Change (e_prune, v_prune, tree);
735	}
736	}
737	}
738
739	/* If the best candidate has a positive estimated change in
740	** likelihood, perform that move.
741	*/
742	if (optim_cand[0]->delta_lk > 1.0/BIG)
743	{
744	improvement = 1;
745	Make_Move (optim_cand[0],0,tree);
746	}
747
748	/*
749	** If there was no improvement at all, try local t_edge length optimization at the
750	** regraft position.
751	*/
752	if (!improvement)
753	{
754	/* PhyML_Printf ("\n. - performing local t_edge length optimizations"); */
755	if ((candidate = Find_Optim_Local (tree)) >= 0)
756	{
757	improvement = 1;
758	Make_Move (optim_cand[candidate],1,tree);
759	}
760	}
761
762	/*
763	** If there was still no improvement, try global t_edge length optimization.
764	*/
765	if (!improvement)
766	{
767	/* PhyML_Printf ("\n. - performing global t_edge length optimization"); */
768	if ((candidate = Find_Optim_Globl (tree)) >= 0)
769	{
770	improvement = 1;
771	Make_Move (optim_cand[candidate],2,tree);
772	}
773	}
774
775	/*
776	** Optimize all t_edge lengths again to make sure we got an updated
777	** likelihood value.
778	*/
779	Set_Both_Sides(YES,tree);
780	cur_lk = Lk(NULL,tree);
781	root = tree->a_nodes[0];
782	Optimize_Br_Len_Serie (tree);
783	Set_Both_Sides(YES,tree);
784	cur_lk = Lk(NULL,tree);
785	time(&(tree->t_current));
786	if(tree->mod->s_opt->print) Print_Lk(tree,"topology");
787
788	/*
789	** Return the result.
790	*/
791	return (improvement);
792	}
793
794
795	/*
796	** Perform_One_Moves: Perform a round of SPR moves on the tree. Prune each subtree in
797	** turn and calculate the change in tree length for each candidate
798	** regraft position. Estimate change in likelihood for the most
799	** promising moves, and perform the first move that results in an
800	** improvement. If no improvements were found at all, try local edge
801	** length optimization. If still no improvement, try global edge
802	** length optimization.
803	**
804	** Parameters:
805	** - tree: The tree to perform the SPR moves on.
806	** - max_size: The maximum size (= number of taxa) of the subtrees to be
807	** pruned. If m=0 or m>ntax, all possible prunings will be
808	** considered.
809	**
810	** Returns:
811	** If the current tree could be improved: 1.
812	** Otherwise: 0.
813	*/
814
815	int Perform_One_SPR(t_tree *tree, int max_size)
816	{
817	int nr_edges, i, j, candidate, improvement;
818	t_node root, v_prune;
819	t_edge *e_prune;
820
821	/*
822	** Calculate the average subtree distances.
823	*/
824
825	root = tree->a_nodes[0];
826	PostOrder_v (tree, root->v[2], root->b[2]);
827
828	/*
829	** Initialize the array of optimization candidates.
830	*/
831	for (i = 0; i < tree->mod->s_opt->wim_n_optim; i++)
832	{
833	optim_cand[i]->delta_lk = -1.0*BIG;
834	optim_cand[i]->d_L = -1.0*BIG;
835	}
836
837	/*
838	** Try all possible SPR moves and perform the ones that give an improvement.
839	*/
840	nr_edges = 2*tree->n_otu - 3;
841	cur_lk = tree->c_lnL;
842	improvement = 0;
843
844	/* PhyML_Printf("\n >>>>>>>>>>>>>>>>>>"); */
845	/* PhyML_Printf("\n. cur_lk = %f %f",cur_lk,Lk(tree)); */
846
847	/* PhyML_Printf ("\n. Trying SPR moves"); */
848	/* PhyML_Printf ("\n. - calculating tree distances and estimating likelihoods"); */
849
850	for(i = 0; i < nr_edges; i++)
851	{
852	/*
853	** Get the next prune edge.
854	*/
855	e_prune = tree->a_edges[i];
856	/*
857	** Try right subtree if appropriate.
858	*/
859	if (!e_prune->left->tax)
860	{
861	/*
862	** Clear the regraft candidate list.
863	*/
864	for (j = 0; j < tree->mod->s_opt->wim_n_rgrft; j++)
865	{
866	rgrft_cand[j]->d_L = -1.0*BIG;
867	}
868	v_prune = e_prune->left;
869	/* if ((max_size == 0) \|\| (e_prune->num_tax_rght <= max_size)) */
870	{
871	/*
872	** Calculate changes in tree length, and estimate changes in likelihood for
873	** the most promising candidates. Perform moves that give an improvement.
874	*/
875	Calc_Tree_Length (e_prune, v_prune, tree);
876	if ((candidate = Est_Lk_Change (e_prune, v_prune, tree)) >= 0)
877	{
878	improvement = 1;
879	Make_Move (rgrft_cand[candidate],0,tree);
880	}
881	}
882	}
883	/*
884	** Try left subtree if appropriate.
885	*/
886	if (!e_prune->rght->tax && !improvement)
887	{
888	/*
889	** Clear the regraft candidate list.
890	*/
891	for (j = 0; j < tree->mod->s_opt->wim_n_rgrft; j++)
892	{
893	rgrft_cand[j]->d_L = -1.0*BIG;
894	}
895	v_prune = e_prune->rght;
896	/* if ((max_size == 0) \|\| (e_prune->num_tax_left <= max_size)) */
897	{
898	/*
899	** Calculate changes in tree length, and estimate changes in likelihood for
900	** the most promising candidates. Perform moves that give an improvement.
901	*/
902	Calc_Tree_Length (e_prune, v_prune, tree);
903	if ((candidate = Est_Lk_Change (e_prune, v_prune, tree)) >= 0)
904	{
905	improvement = 1;
906	Make_Move (rgrft_cand[candidate],0,tree);
907	/* PhyML_Printf("\n. Make simple move"); */
908	/* PhyML_Printf("\n. lk after simple move = %f",Lk(tree)); */
909	}
910	}
911	}
912	if(improvement) break;
913	}
914
915	/*
916	** If there was no improvement at all, try local t_edge length optimization at the
917	** regraft position.
918	*/
919
920	/* PhyML_Printf("\n. before local = %f %f",tree->c_lnL,Lk(tree)); */
921
922	if (!improvement)
923	{
924	/* PhyML_Printf ("\n. - performing local t_edge length optimizations"); */
925	if ((candidate = Find_Optim_Local (tree)) >= 0)
926	{
927	/* PhyML_Printf("\n. make local move"); */
928	improvement = 1;
929	Make_Move (optim_cand[candidate],1,tree);
930	/* PhyML_Printf("\n. lk after local move = %f",Lk(tree)); */
931	}
932	}
933
934
935	/*
936	** If there was still no improvement, try global t_edge length optimization.
937	*/
938	/* PhyML_Printf("\n. before global = %f %f",tree->c_lnL,Lk(tree)); */
939
940	if (!improvement)
941	{
942	/* PhyML_Printf ("\n. - performing global t_edge length optimization"); */
943	if ((candidate = Find_Optim_Globl (tree)) >= 0)
944	{
945	/* PhyML_Printf("\n. make global move"); */
946	improvement = 1;
947	Make_Move (optim_cand[candidate],2,tree);
948	/* PhyML_Printf("\n. lk after global move = %f",Lk(tree)); */
949	}
950	}
951
952	/* PhyML_Printf("\n. after all = %f %f",tree->c_lnL,Lk(tree)); */
953
954	/*
955	** Optimize all t_edge lengths again to make sure we got an updated
956	** likelihood value.
957	*/
958	Set_Both_Sides(YES,tree);
959	cur_lk = Lk(NULL,tree);
960	root = tree->a_nodes[0];
961	Optimize_Br_Len_Serie (tree);
962	Set_Both_Sides(YES,tree);
963	cur_lk = Lk(NULL,tree);
964	time(&(tree->t_current));
965	if(tree->mod->s_opt->print) Print_Lk(tree,"topology");
966
967	/*
968	** Return the result.
969	*/
970	return (improvement);
971	}
972
973
974	/*
975	** Calc_Tree_Length: Calculate the change in tree length, given a pruned subtree,
976	** for each possible regraft position.
977	**
978	** Parameters:
979	** - e_prune: The t_edge at which the subtree is pruned.
980	** - v_prune: The root of the pruned subtree.
981	*/
982
983	void Calc_Tree_Length (t_edge e_prune, t_node v_prune, t_tree *tree)
984	{
985	int i, d0, d1, d2;
986	phydbl d_uu;
987	t_node u_prune, u1, *u2;
988
989	/*
990	** Get the directions from t_node v_prune.
991	*/
992	d0 = -1;
993	u_prune = NULL;
994	for (i = 0; i < 3; i++)
995	{
996	if (v_prune->b[i] == e_prune)
997	{
998	d0 = i;
999	u_prune = v_prune->v[i];
1000	break;
1001	}
1002	}
1003	d1 = (d0 + 1) % 3;
1004	d2 = 3 - d0 - d1;
1005
1006	/*
1007	** Get the relevant average subtree distance within the pruned subtree.
1008	*/
1009	if (!u_prune->tax)
1010	{
1011	u1 = u2 = NULL;
1012	for (i = 0; i < 3; i++)
1013	{
1014	if (u_prune->b[i] != e_prune)
1015	{
1016	if (u1 == NULL)
1017	{
1018	u1 = u_prune->v[i];
1019	}
1020	else
1021	{
1022	u2 = u_prune->v[i];
1023	}
1024	}
1025	}
1026	d_uu = subtree_dist[u1->num][u2->num];
1027	}
1028	else
1029	{
1030	d_uu = 0.0;
1031	}
1032
1033	/*
1034	** Recursively calculate the change in tree length for each
1035	** possible regraft position.
1036	**
1037	** First recurse into direction d1.
1038	*/
1039	if (!v_prune->v[d1]->tax)
1040	{
1041	u1 = u2 = NULL;
1042	for (i = 0; i < 3; i++)
1043	{
1044	if (v_prune->v[d1]->b[i] != v_prune->b[d1])
1045	{
1046	if (u1 == NULL)
1047	{
1048	u1 = v_prune->v[d1]->v[i];
1049	}
1050	else
1051	{
1052	u2 = v_prune->v[d1]->v[i];
1053	}
1054	}
1055	}
1056	Tree_Length(v_prune, u_prune, v_prune->v[d1], v_prune->v[d2], u1, v_prune->v[d2],
1057	u2, subtree_dist[u_prune->num][v_prune->v[d2]->num], d_uu, 0.0, 1, tree);
1058	Tree_Length(v_prune, u_prune, v_prune->v[d1], v_prune->v[d2], u2, v_prune->v[d2],
1059	u1, subtree_dist[u_prune->num][v_prune->v[d2]->num], d_uu, 0.0, 1, tree);
1060	}
1061	/*
1062	** Next recurse into direction d2.
1063	*/
1064	if (!v_prune->v[d2]->tax)
1065	{
1066	u1 = u2 = NULL;
1067	for (i = 0; i < 3; i++)
1068	{
1069	if (v_prune->v[d2]->b[i] != v_prune->b[d2])
1070	{
1071	if (u1 == NULL)
1072	{
1073	u1 = v_prune->v[d2]->v[i];
1074	}
1075	else
1076	{
1077	u2 = v_prune->v[d2]->v[i];
1078	}
1079	}
1080	}
1081	Tree_Length(v_prune, u_prune, v_prune->v[d2], v_prune->v[d1], u1, v_prune->v[d1],
1082	u2, subtree_dist[u_prune->num][v_prune->v[d1]->num], d_uu, 0.0, 1, tree);
1083	Tree_Length(v_prune, u_prune, v_prune->v[d2], v_prune->v[d1], u2, v_prune->v[d1],
1084	u1, subtree_dist[u_prune->num][v_prune->v[d1]->num], d_uu, 0.0, 1, tree);
1085	}
1086	}
1087
1088
1089	/*
1090	** Tree_Length: Recursively calculate the change in tree length for a given pruned
1091	** subtree and regraft position.
1092	**
1093	** Parameters:
1094	** - v_prune: The root of the pruned subtree.
1095	** - u_prune: The t_node adjacent to v_p along the pruned edge.
1096	** - v_n: The t_node adjacent to the regraft t_edge in the "backward" direction.
1097	** - v_n_1: The previous v_n.
1098	** - v_nx1: The t_node adjacent to the regrafting t_edge in the "forward" direction.
1099	** - v_0: The other t_node originally adjacent to v_p;
1100	** - u_n: The other t_node adjecent to v_n (besides v_n_1 and v_nx1);
1101	** - d_uv_1: The distance between u_p and v_n_1;
1102	** - d_uu: The subtree distance between descendants of u_prune.
1103	** - d_L_1: The previous change in tree length.
1104	** - n: The current distance from the prune position.
1105	*/
1106
1107	void Tree_Length (t_node v_prune, t_node u_prune, t_node v_n, t_node v_n_1, t_node *v_nx1,
1108	t_node v_0, t_node u_n, phydbl d_up_v_1, phydbl d_L_1, phydbl d_uu,
1109	int n, t_tree *tree)
1110	{
1111	int i, j;
1112	phydbl d_un_v, d_up_v, d_L;
1113	t_node u1, u2;
1114	t_edge e_prune, e_regraft;
1115	_move_ *tmp_cand;
1116
1117	/*
1118	** Update the path and number of calculations.
1119	*/
1120	path[n] = v_n;
1121	nr_d_L++;
1122	e_prune = NULL;
1123	e_regraft = NULL;
1124
1125	/*
1126	** Calculate the change in tree length for the current pruned subtree and regraft
1127	** position.
1128	*/
1129	if (n == 1)
1130	{
1131	d_un_v = subtree_dist[u_n->num][v_0->num];
1132	}
1133	else
1134	{
1135	d_un_v = subtree_dist[u_n->num][v_n_1->num] -
1136	(pow (0.5, n) * subtree_dist[u_n->num][u_prune->num]) +
1137	(pow (0.5, n) * subtree_dist[u_n->num][v_0->num]);
1138	}
1139	d_up_v = 0.5 * (d_up_v_1 + subtree_dist[u_prune->num][u_n->num]);
1140	/*
1141	** Alternative method for calculating d_up_v. Just kept it around for reference...
1142	**
1143	d_up_v = subtree_dist[u_prune->num][u_n->num] - (0.5 * d_uu);
1144	if (!u_n->tax)
1145	{
1146	u1 = u2 = NULL;
1147	for (i = 0; i < 3; i++)
1148	{
1149	if (u_n->v[i] != v_n)
1150	{
1151	if (u1 == NULL)
1152	{
1153	u1 = u_n->v[i];
1154	}
1155	else
1156	{
1157	u2 = u_n->v[i];
1158	}
1159	}
1160	}
1161	d_up_v -= 0.5 * subtree_dist[u1->num][u2->num];
1162	}
1163	for (i = 0; i < 3; i++)
1164	{
1165	if (u_n->v[i] == v_n)
1166	{
1167	d_up_v -= u_n->b[i]->l->v;
1168	break;
1169	}
1170	}
1171	*/
1172	d_L = d_L_1 + 0.25*((d_up_v_1 + subtree_dist[u_n->num][v_nx1->num]) -
1173	(d_un_v + subtree_dist[u_prune->num][v_nx1->num]));
1174
1175	/*
1176	** If the change is within the tree->mod->s_opt->wim_n_rgrft best ones so far, save it.
1177	*/
1178	if (d_L > rgrft_cand[tree->mod->s_opt->wim_n_rgrft-1]->d_L)
1179	{
1180	for (i = 0; i < 3; i++)
1181	{
1182	if (v_prune->v[i] == u_prune)
1183	{
1184	e_prune = v_prune->b[i];
1185	}
1186	if (v_n->v[i] == v_nx1)
1187	{
1188	e_regraft = v_n->b[i];
1189	}
1190	}
1191	i = tree->mod->s_opt->wim_n_rgrft-1;
1192	rgrft_cand[i]->v_prune = v_prune;
1193	rgrft_cand[i]->u_prune = u_prune;
1194	rgrft_cand[i]->v_n = v_n;
1195	rgrft_cand[i]->v_nx1 = v_nx1;
1196	rgrft_cand[i]->u_n = u_n;
1197	rgrft_cand[i]->e_prune = e_prune;
1198	rgrft_cand[i]->e_regraft = e_regraft;
1199	rgrft_cand[i]->d_L = d_L;
1200	rgrft_cand[i]->d_up_v = d_up_v;
1201	rgrft_cand[i]->d_un_v = d_un_v;
1202	rgrft_cand[i]->dist = n;
1203	for (j = 1; j <= n; j++)
1204	{
1205	rgrft_cand[i]->path[j] = path[j];
1206	}
1207
1208	rgrft_cand[i]->path[n+1] = v_nx1;
1209	/*
1210	** Move the candidate to the appropriate position in the list, so the list
1211	** remains sorted in decreasing d_L value.
1212	*/
1213	while ((i > 0) && (rgrft_cand[i]->d_L > rgrft_cand[i-1]->d_L))
1214	{
1215	tmp_cand = rgrft_cand[i];
1216	rgrft_cand[i] = rgrft_cand[i-1];
1217	rgrft_cand[i-1] = tmp_cand;
1218	i--;
1219	}
1220	}
1221
1222	/*
1223	** Recurse.
1224	*/
1225	if (n < tree->mod->s_opt->wim_max_dist)
1226	{
1227	if (!v_nx1->tax)
1228	{
1229	u1 = u2 = NULL;
1230	for (i = 0; i < 3; i++)
1231	{
1232	if (v_nx1->v[i] != v_n)
1233	{
1234	if (u1 == NULL)
1235	{
1236	u1 = v_nx1->v[i];
1237	}
1238	else
1239	{
1240	u2 = v_nx1->v[i];
1241	}
1242	}
1243	}
1244	Tree_Length (v_prune, u_prune, v_nx1, v_n, u1, v_0, u2, d_up_v, d_uu, d_L, n+1, tree);
1245	Tree_Length (v_prune, u_prune, v_nx1, v_n, u2, v_0, u1, d_up_v, d_uu, d_L, n+1, tree);
1246	}
1247	}
1248	}
1249
1250
1251	/*
1252	** Est_Lk_Change: Estimate the changes in likelihood for the most promising candidate
1253	** regraft positions given a pruned subtree.
1254	**
1255	** Parameters:
1256	** - e_prune: The t_edge at which the subtree was pruned.
1257	** - v_prune: The root of the pruned subtree.
1258	** - tree: The tree on which to do the calculations.
1259	**
1260	** Returns:
1261	** If an improvement as found: The candidate which gives the improvement (which
1262	** will be the first one found).
1263	** Otherwise: -1.
1264	*/
1265
1266	int Est_Lk_Change (t_edge e_prune, t_node v_prune, t_tree *tree)
1267	{
1268	int i, j, cand, best_cand, d0, d1, d2, n, pat, cat, ste;
1269	phydbl d_uu, best_d_lk, l_connect, l_01, l_02, l_12, l_est[3], new_lk,
1270	l_simple[3], l_dist[3];
1271	phydbl p_lk1_tmp, p_lk2_tmp, *p_lk;
1272	int *p_sum;
1273	t_node u_prune, v_n, v_nx1, u1, *u2;
1274	t_edge e_regraft, e_tmp;
1275	_move_ *tmp_cand;
1276	int dim1, dim2;
1277
1278
1279	dim1 = tree->mod->ns * tree->mod->ras->n_catg;
1280	dim2 = tree->mod->ras->n_catg;
1281
1282	/*
1283	** Get the directions from t_node v_prune.
1284	*/
1285	d0 = -1;
1286	u_prune = NULL;
1287	for (i = 0; i < 3; i++)
1288	{
1289	if (v_prune->b[i] == e_prune)
1290	{
1291	d0 = i;
1292	u_prune = v_prune->v[i];
1293	break;
1294	}
1295	}
1296	d1 = (d0 + 1) % 3;
1297	d2 = 3 - d0 - d1;
1298
1299	/*
1300	** Copy the relevant partial likelihoods to the temporary regraft structure.
1301	** We can point to the original matrices, cos they won't be changed anyway.
1302	*/
1303	if (v_prune == e_prune->left)
1304	{
1305	v_tmp->b[0]->p_lk_rght = e_prune->p_lk_rght;
1306	v_tmp->b[0]->sum_scale_rght = e_prune->sum_scale_rght;
1307	}
1308	else
1309	{
1310	v_tmp->b[0]->p_lk_rght = e_prune->p_lk_left;
1311	v_tmp->b[0]->sum_scale_rght = e_prune->sum_scale_left;
1312	}
1313	v_tmp->num = v_prune->num;
1314	v_tmp->v[0]->num = u_prune->num;
1315	v_tmp->b[0]->num = e_prune->num;
1316
1317	/*
1318	** Estimate the length of the t_edge that will connect the two "detached" nodes
1319	** after pruning. (The average of the sum of the lengths of the original two
1320	** edges and the average subtree distance based estimate.)
1321	*/
1322	l_connect = subtree_dist[v_prune->v[d1]->num][v_prune->v[d2]->num];
1323	if (!v_prune->v[d1]->tax)
1324	{
1325	u1 = u2 = NULL;
1326	for (i = 0; i < 3; i++)
1327	{
1328	if (v_prune->v[d1]->b[i] != v_prune->b[d1])
1329	{
1330	if (u1 == NULL)
1331	{
1332	u1 = v_prune->v[d1]->v[i];
1333	}
1334	else
1335	{
1336	u2 = v_prune->v[d1]->v[i];
1337	}
1338	}
1339	}
1340	l_connect -= 0.5 * subtree_dist[u1->num][u2->num];
1341	}
1342	if (!v_prune->v[d2]->tax)
1343	{
1344	u1 = u2 = NULL;
1345	for (i = 0; i < 3; i++)
1346	{
1347	if (v_prune->v[d2]->b[i] != v_prune->b[d2])
1348	{
1349	if (u1 == NULL)
1350	{
1351	u1 = v_prune->v[d2]->v[i];
1352	}
1353	else
1354	{
1355	u2 = v_prune->v[d2]->v[i];
1356	}
1357	}
1358	}
1359	l_connect -= 0.5 * subtree_dist[u1->num][u2->num];
1360	}
1361	l_connect += (v_prune->b[d1]->l->v + v_prune->b[d2]->l->v);
1362	l_connect /= 2.0;
1363
1364	/*
1365	** Temporarily swap the relevant partial likelihoods at the prune site.
1366	**
1367	** Direction d1.
1368	*/
1369	if (v_prune == v_prune->b[d1]->left)
1370	{
1371	p_lk1_tmp = v_prune->b[d1]->p_lk_left;
1372	if (v_prune == v_prune->b[d2]->left)
1373	{
1374	v_prune->b[d1]->p_lk_left = v_prune->b[d2]->p_lk_rght;
1375	}
1376	else
1377	{
1378	v_prune->b[d1]->p_lk_left = v_prune->b[d2]->p_lk_left;
1379	}
1380	}
1381	else
1382	{
1383	p_lk1_tmp = v_prune->b[d1]->p_lk_rght;
1384	if (v_prune == v_prune->b[d2]->left)
1385	{
1386	v_prune->b[d1]->p_lk_rght = v_prune->b[d2]->p_lk_rght;
1387	}
1388	else
1389	{
1390	v_prune->b[d1]->p_lk_rght = v_prune->b[d2]->p_lk_left;
1391	}
1392	}
1393	/*
1394	** Direction d2.
1395	*/
1396	if (v_prune == v_prune->b[d2]->left)
1397	{
1398	p_lk2_tmp = v_prune->b[d2]->p_lk_left;
1399	if (v_prune == v_prune->b[d1]->left)
1400	{
1401	v_prune->b[d2]->p_lk_left = v_prune->b[d1]->p_lk_rght;
1402	}
1403	else
1404	{
1405	v_prune->b[d2]->p_lk_left = v_prune->b[d1]->p_lk_left;
1406	}
1407	}
1408	else
1409	{
1410	p_lk2_tmp = v_prune->b[d2]->p_lk_rght;
1411	if (v_prune == v_prune->b[d1]->left)
1412	{
1413	v_prune->b[d2]->p_lk_rght = v_prune->b[d1]->p_lk_rght;
1414	}
1415	else
1416	{
1417	v_prune->b[d2]->p_lk_rght = v_prune->b[d1]->p_lk_left;
1418	}
1419	}
1420
1421	/*
1422	** Temporarily set the t_edge lengths and update transition prob's at the
1423	** prune site.
1424	*/
1425	v_prune->b[d1]->l_old->v = v_prune->b[d1]->l->v;
1426	v_prune->b[d2]->l_old->v = v_prune->b[d2]->l->v;
1427	v_prune->b[d1]->l->v = l_connect;
1428	v_prune->b[d2]->l->v = l_connect;
1429	Update_PMat_At_Given_Edge (v_prune->b[d1], tree);
1430	Update_PMat_At_Given_Edge (v_prune->b[d2], tree);
1431
1432	/*
1433	** Get the relevant average subtree distance within the pruned subtree.
1434	*/
1435	if (!u_prune->tax)
1436	{
1437	u1 = u2 = NULL;
1438	for (i = 0; i < 3; i++)
1439	{
1440	if (u_prune->b[i] != e_prune)
1441	{
1442	if (u1 == NULL)
1443	{
1444	u1 = u_prune->v[i];
1445	}
1446	else
1447	{
1448	u2 = u_prune->v[i];
1449	}
1450	}
1451	}
1452	d_uu = subtree_dist[u1->num][u2->num];
1453	}
1454	else
1455	{
1456	d_uu = 0.0;
1457	}
1458
1459	/*
1460	** Try each candidate SPR and estimate the change in likelihood.
1461	*/
1462	best_d_lk = 1.0/BIG;
1463	best_cand = -1;
1464	for (cand = 0; cand < tree->mod->s_opt->wim_n_rgrft; cand++)
1465	{
1466	/*
1467	** If there are no more candidates, bail out...
1468	*/
1469	if (FABS(rgrft_cand[cand]->d_L - 1.0*BIG) < SMALL)
1470	{
1471	break;
1472	}
1473	else
1474	{
1475	nr_d_lk++;
1476	}
1477
1478	/*
1479	** Get the relevant nodes and edges.
1480	*/
1481	v_n = rgrft_cand[cand]->v_n;
1482	v_nx1 = rgrft_cand[cand]->v_nx1;
1483	e_regraft = rgrft_cand[cand]->e_regraft;
1484
1485	/*
1486	** Update the relevant partial likelihoods along the path between the prune
1487	** and regraft positions (temporarily save the first one).
1488	*/
1489	n = rgrft_cand[cand]->dist;
1490	e_tmp = NULL;
1491	p_lk = NULL;
1492	p_sum = NULL;
1493	for (i = 1; i <= n; i++)
1494	{
1495	/*
1496	** Get the next t_edge along the path.
1497	*/
1498	for (j = 0; j < 3; j++)
1499	{
1500	if (rgrft_cand[cand]->path[i]->v[j] == rgrft_cand[cand]->path[i+1])
1501	{
1502	e_tmp = rgrft_cand[cand]->path[i]->b[j];
1503	break;
1504	}
1505	}
1506	if (i == 1)
1507	{
1508	/*
1509	** Save the first partial likelihood along the path.
1510	*/
1511	if (rgrft_cand[cand]->path[i] == e_tmp->left)
1512	{
1513	p_lk = e_tmp->p_lk_left;
1514	p_sum = e_tmp->sum_scale_left;
1515	}
1516	else
1517	{
1518	p_lk = e_tmp->p_lk_rght;
1519	p_sum = e_tmp->sum_scale_rght;
1520	}
1521
1522	for (pat = 0; pat < tree->n_pattern; pat++)
1523	{
1524	sum_scale_tmp[pat] = p_sum[pat];
1525	for (cat = 0; cat < tree->mod->ras->n_catg; cat++)
1526	{
1527	for (ste = 0; ste < tree->mod->ns; ste++)
1528	{
1529	p_lk_tmp[patdim1+catdim2+ste] = p_lk[patdim1+catdim2+ste];
1530	}
1531	}
1532	}
1533	}
1534	Update_P_Lk (tree, e_tmp, rgrft_cand[cand]->path[i]);
1535	}
1536	if (v_n == e_regraft->left)
1537	{
1538	v_tmp->b[1]->p_lk_rght = e_regraft->p_lk_left;
1539	v_tmp->b[2]->p_lk_rght = e_regraft->p_lk_rght;
1540	v_tmp->b[1]->sum_scale_rght = e_regraft->sum_scale_left;
1541	v_tmp->b[2]->sum_scale_rght = e_regraft->sum_scale_rght;
1542	}
1543	else
1544	{
1545	v_tmp->b[1]->p_lk_rght = e_regraft->p_lk_rght;
1546	v_tmp->b[2]->p_lk_rght = e_regraft->p_lk_left;
1547	v_tmp->b[1]->sum_scale_rght = e_regraft->sum_scale_rght;
1548	v_tmp->b[2]->sum_scale_rght = e_regraft->sum_scale_left;
1549	}
1550
1551	/*
1552	** Estimate t_edge lengths of the three relevant regraft edges based on
1553	** average subtree distances.
1554	**
1555	** l_01
1556	*/
1557	/*
1558	** Alternative method of estimating l_01. Kept it around for reference...
1559	**
1560	l_01 = subtree_dist[u_prune->num][u_n->num] - (0.5 * d_uu);
1561	if (!u_n->tax)
1562	{
1563	u1 = u2 = NULL;
1564	for (i = 0; i < 3; i++)
1565	{
1566	if (u_n->v[i] != v_n)
1567	{
1568	if (u1 == NULL)
1569	{
1570	u1 = u_n->v[i];
1571	}
1572	else
1573	{
1574	u2 = u_n->v[i];
1575	}
1576	}
1577	}
1578	l_01 -= 0.5 * subtree_dist[u1->num][u2->num];
1579	}
1580	for (i = 0; i < 3; i++)
1581	{
1582	if (u_n->v[i] == v_n)
1583	{
1584	l_01 -= u_n->b[i]->l->v;
1585	break;
1586	}
1587	}
1588	*/
1589	l_01 = rgrft_cand[cand]->d_up_v - (0.5 * rgrft_cand[cand]->d_un_v) -
1590	(0.5 * d_uu);
1591	/*
1592	** l_02
1593	*/
1594	l_02 = subtree_dist[u_prune->num][v_nx1->num] - (0.5 * d_uu);
1595	if (!v_nx1->tax)
1596	{
1597	u1 = u2 = NULL;
1598	for (i = 0; i < 3; i++)
1599	{
1600	if (v_nx1->v[i] != v_n)
1601	{
1602	if (u1 == NULL)
1603	{
1604	u1 = v_nx1->v[i];
1605	}
1606	else
1607	{
1608	u2 = v_nx1->v[i];
1609	}
1610	}
1611	}
1612	l_02 -= (0.5 * subtree_dist[u1->num][u2->num]);
1613	}
1614	/*
1615	** l_12
1616	*/
1617	l_12 = e_regraft->l->v;
1618	/*
1619	** Simple estimates.
1620	*/
1621	l_simple[0] = l_02 - (0.5*e_regraft->l->v);
1622	l_simple[1] = 0.5 * e_regraft->l->v;
1623	l_simple[2] = 0.5 * e_regraft->l->v;
1624	/*
1625	** Average subtree distance based estimates.
1626	*/
1627	l_dist[0] = 0.5 * ( l_01 + l_02 - l_12);
1628	l_dist[1] = 0.5 * ( l_01 - l_02 + l_12);
1629	l_dist[2] = 0.5 * (-l_01 + l_02 + l_12);
1630	/*
1631	** Take the average of the two estimates.
1632	*/
1633	l_est[0] = (l_simple[0] + l_dist[0]) / 2.0;
1634	l_est[1] = (l_simple[1] + l_dist[1]) / 2.0;
1635	l_est[2] = (l_simple[2] + l_dist[2]) / 2.0;
1636
1637	/*
1638	** Set the t_edge lengths and update the relevant transition prob's and
1639	** partial likelihoods in the temporary regraft structure.
1640	*/
1641
1642	for (i = 0; i < 3; i++)
1643	{
1644	v_tmp->b[i]->l->v = l_est[i]; /* TO DO */
1645	Update_PMat_At_Given_Edge (v_tmp->b[i], tree);
1646	}
1647
1648	/* Beg SG 18 May 2007 */
1649	if(tree->mod->s_opt->wim_inside_opt)
1650	{
1651	Triple_Dist(v_tmp,tree,0);
1652	For(i,3) l_est[i] = v_tmp->b[i]->l->v;
1653	}
1654	/* End SG 18 May 2007 */
1655
1656
1657	/*
1658	** Calculate the change in likelihood locally. Save it and the estimated edge
1659	** lengths in the current candidate in the list.
1660	*/
1661	Update_P_Lk (tree, v_tmp->b[0], v_tmp);
1662	new_lk = Lk(v_tmp->b[0],tree);
1663	/* PhyML_Printf("\n. new_lk = %f",new_lk); */
1664
1665	rgrft_cand[cand]->delta_lk = new_lk - cur_lk;
1666	rgrft_cand[cand]->rgrft_rank = cand;
1667	rgrft_cand[cand]->optim_rank = -1;
1668	rgrft_cand[cand]->globl_rank = -1;
1669	rgrft_cand[cand]->l_connect = l_connect;
1670	for (i = 0; i < 3; i++)
1671	{
1672	rgrft_cand[cand]->l_est[i] = v_tmp->b[i]->l->v;
1673	}
1674	if (rgrft_cand[cand]->delta_lk > best_d_lk)
1675	{
1676	best_d_lk = rgrft_cand[cand]->delta_lk;
1677	best_cand = cand;
1678	}
1679
1680	/*
1681	** If the change is within the tree->mod->s_opt->wim_n_optim best ones, save it in the list of
1682	** optimization candidates.
1683	*/
1684	if (rgrft_cand[cand]->delta_lk > optim_cand[tree->mod->s_opt->wim_n_optim-1]->delta_lk)
1685	{
1686	i = tree->mod->s_opt->wim_n_optim-1;
1687	optim_cand[i]->v_prune = rgrft_cand[cand]->v_prune;
1688	optim_cand[i]->u_prune = rgrft_cand[cand]->u_prune;
1689	optim_cand[i]->v_n = rgrft_cand[cand]->v_n;
1690	optim_cand[i]->v_nx1 = rgrft_cand[cand]->v_nx1;
1691	optim_cand[i]->u_n = rgrft_cand[cand]->u_n;
1692	optim_cand[i]->e_prune = rgrft_cand[cand]->e_prune;
1693	optim_cand[i]->e_regraft = rgrft_cand[cand]->e_regraft;
1694	optim_cand[i]->d_L = rgrft_cand[cand]->d_L;
1695	optim_cand[i]->dist = rgrft_cand[cand]->dist;
1696	optim_cand[i]->rgrft_rank = rgrft_cand[cand]->rgrft_rank;
1697	optim_cand[i]->optim_rank = rgrft_cand[cand]->optim_rank;
1698	optim_cand[i]->globl_rank = rgrft_cand[cand]->globl_rank;
1699	optim_cand[i]->l_connect = rgrft_cand[cand]->l_connect;
1700	for (j = 0; j < 3; j++)
1701	{
1702	optim_cand[i]->l_est[j] = rgrft_cand[cand]->l_est[j];
1703	}
1704	optim_cand[i]->delta_lk = rgrft_cand[cand]->delta_lk;
1705	/*
1706	** Move the candidate to the appropriate position in the list, so the list
1707	** remains sorted in decreasing delta_Lk value.
1708	*/
1709	while ((i > 0) && (optim_cand[i]->delta_lk > optim_cand[i-1]->delta_lk))
1710	{
1711	tmp_cand = optim_cand[i];
1712	optim_cand[i] = optim_cand[i-1];
1713	optim_cand[i-1] = tmp_cand;
1714	i--;
1715	}
1716	}
1717
1718	/*
1719	** Reset the partial likelihoods along the path.
1720	*/
1721	for (pat = 0; pat < tree->n_pattern; pat++)
1722	{
1723	p_sum[pat] = sum_scale_tmp[pat];
1724	for (cat = 0; cat < tree->mod->ras->n_catg; cat++)
1725	{
1726	for (ste = 0; ste < tree->mod->ns; ste++)
1727	{
1728	p_lk[patdim1+catdim2+ste] = p_lk_tmp[patdim1+catdim2+ste];
1729	}
1730	}
1731	}
1732	n = rgrft_cand[cand]->dist;
1733	for (i = 2; i <= n; i++)
1734	{
1735	for (j = 0; j < 3; j++)
1736	{
1737	if (rgrft_cand[cand]->path[i]->v[j] == rgrft_cand[cand]->path[i+1])
1738	{
1739	e_tmp = rgrft_cand[cand]->path[i]->b[j];
1740	break;
1741	}
1742	}
1743	Update_P_Lk (tree, e_tmp, rgrft_cand[cand]->path[i]);
1744	}
1745
1746	/*
1747	** If an improvement was found, forget the other candidates...
1748	*/
1749	if (best_cand >= 0)
1750	{
1751	break;
1752	}
1753	}
1754
1755	/*
1756	** Swap back the relevant partial likelihoods at the prune site.
1757	*/
1758	if (v_prune == v_prune->b[d1]->left)
1759	{
1760	v_prune->b[d1]->p_lk_left = p_lk1_tmp;
1761	}
1762	else
1763	{
1764	v_prune->b[d1]->p_lk_rght = p_lk1_tmp;
1765	}
1766	if (v_prune == v_prune->b[d2]->left)
1767	{
1768	v_prune->b[d2]->p_lk_left = p_lk2_tmp;
1769	}
1770	else
1771	{
1772	v_prune->b[d2]->p_lk_rght = p_lk2_tmp;
1773	}
1774
1775	/*
1776	** Reset the relevant t_edge lengths and transition prob's at the prune site.
1777	*/
1778	v_prune->b[d1]->l->v = v_prune->b[d1]->l_old->v;
1779	v_prune->b[d2]->l->v = v_prune->b[d2]->l_old->v;
1780	Update_PMat_At_Given_Edge (v_prune->b[d1], tree);
1781	Update_PMat_At_Given_Edge (v_prune->b[d2], tree);
1782
1783	/*
1784	** Return the best candidate.
1785	*/
1786	return (best_cand);
1787	}
1788
1789
1790	/*
1791	** Best_Lk_Change: Estimate the changes in likelihood for the most promising candidate
1792	** regraft positions given a pruned subtree and save the best one.
1793	**
1794	** Parameters:
1795	** - e_prune: The t_edge at which the subtree was pruned.
1796	** - v_prune: The root of the pruned subtree.
1797	** - tree: The tree on which to do the calculations.
1798	**
1799	** Returns:
1800	** The candidate which gives the best (possibly negative) improvement.
1801	*/
1802
1803	int Best_Lk_Change (t_edge e_prune, t_node v_prune, t_tree *tree)
1804	{
1805	int i, j, cand, best_cand, d0, d1, d2, n, pat, cat, ste;
1806	phydbl d_uu, best_d_lk, l_connect, l_01, l_02, l_12, l_est[3], new_lk, l_simple[3], l_dist[3];
1807	phydbl p_lk1_tmp, p_lk2_tmp, *p_lk;
1808	int *p_sum;
1809	t_node u_prune, v_n, v_nx1, u1, *u2;
1810	t_edge e_regraft, e_tmp;
1811	_move_ *tmp_cand;
1812	int dim1, dim2;
1813
1814	dim1 = tree->mod->ns * tree->mod->ras->n_catg;
1815	dim2 = tree->mod->ns ;
1816
1817	/*
1818	** Get the directions from t_node v_prune.
1819	*/
1820	d0 = -1;
1821	u_prune = NULL;
1822	for (i = 0; i < 3; i++)
1823	{
1824	if (v_prune->b[i] == e_prune)
1825	{
1826	d0 = i;
1827	u_prune = v_prune->v[i];
1828	break;
1829	}
1830	}
1831	d1 = (d0 + 1) % 3;
1832	d2 = 3 - d0 - d1;
1833
1834	/*
1835	** Copy the relevant partial likelihoods to the temporary regraft structure.
1836	** We can point to the original matrices, cos they won't be changed anyway.
1837	*/
1838	if (v_prune == e_prune->left)
1839	{
1840	v_tmp->b[0]->p_lk_rght = e_prune->p_lk_rght;
1841	v_tmp->b[0]->sum_scale_rght = e_prune->sum_scale_rght;
1842	}
1843	else
1844	{
1845	v_tmp->b[0]->p_lk_rght = e_prune->p_lk_left;
1846	v_tmp->b[0]->sum_scale_rght = e_prune->sum_scale_left;
1847	}
1848	v_tmp->num = v_prune->num;
1849	v_tmp->v[0]->num = u_prune->num;
1850	v_tmp->b[0]->num = e_prune->num;
1851
1852	/*
1853	** Estimate the length of the t_edge that will connect the two "detached" nodes
1854	** after pruning. (The average of the sum of the lengths of the original two
1855	** edges and the average subtree distance based estimate.)
1856	*/
1857	l_connect = subtree_dist[v_prune->v[d1]->num][v_prune->v[d2]->num];
1858	if (!v_prune->v[d1]->tax)
1859	{
1860	u1 = u2 = NULL;
1861	for (i = 0; i < 3; i++)
1862	{
1863	if (v_prune->v[d1]->b[i] != v_prune->b[d1])
1864	{
1865	if (u1 == NULL)
1866	{
1867	u1 = v_prune->v[d1]->v[i];
1868	}
1869	else
1870	{
1871	u2 = v_prune->v[d1]->v[i];
1872	}
1873	}
1874	}
1875	l_connect -= 0.5 * subtree_dist[u1->num][u2->num];
1876	}
1877	if (!v_prune->v[d2]->tax)
1878	{
1879	u1 = u2 = NULL;
1880	for (i = 0; i < 3; i++)
1881	{
1882	if (v_prune->v[d2]->b[i] != v_prune->b[d2])
1883	{
1884	if (u1 == NULL)
1885	{
1886	u1 = v_prune->v[d2]->v[i];
1887	}
1888	else
1889	{
1890	u2 = v_prune->v[d2]->v[i];
1891	}
1892	}
1893	}
1894	l_connect -= 0.5 * subtree_dist[u1->num][u2->num];
1895	}
1896	l_connect += (v_prune->b[d1]->l->v + v_prune->b[d2]->l->v);
1897	l_connect /= 2.0;
1898
1899	/*
1900	** Temporarily swap the relevant partial likelihoods at the prune site.
1901	**
1902	** Direction d1.
1903	*/
1904	if (v_prune == v_prune->b[d1]->left)
1905	{
1906	p_lk1_tmp = v_prune->b[d1]->p_lk_left;
1907	if (v_prune == v_prune->b[d2]->left)
1908	{
1909	v_prune->b[d1]->p_lk_left = v_prune->b[d2]->p_lk_rght;
1910	}
1911	else
1912	{
1913	v_prune->b[d1]->p_lk_left = v_prune->b[d2]->p_lk_left;
1914	}
1915	}
1916	else
1917	{
1918	p_lk1_tmp = v_prune->b[d1]->p_lk_rght;
1919	if (v_prune == v_prune->b[d2]->left)
1920	{
1921	v_prune->b[d1]->p_lk_rght = v_prune->b[d2]->p_lk_rght;
1922	}
1923	else
1924	{
1925	v_prune->b[d1]->p_lk_rght = v_prune->b[d2]->p_lk_left;
1926	}
1927	}
1928	/*
1929	** Direction d2.
1930	*/
1931	if (v_prune == v_prune->b[d2]->left)
1932	{
1933	p_lk2_tmp = v_prune->b[d2]->p_lk_left;
1934	if (v_prune == v_prune->b[d1]->left)
1935	{
1936	v_prune->b[d2]->p_lk_left = v_prune->b[d1]->p_lk_rght;
1937	}
1938	else
1939	{
1940	v_prune->b[d2]->p_lk_left = v_prune->b[d1]->p_lk_left;
1941	}
1942	}
1943	else
1944	{
1945	p_lk2_tmp = v_prune->b[d2]->p_lk_rght;
1946	if (v_prune == v_prune->b[d1]->left)
1947	{
1948	v_prune->b[d2]->p_lk_rght = v_prune->b[d1]->p_lk_rght;
1949	}
1950	else
1951	{
1952	v_prune->b[d2]->p_lk_rght = v_prune->b[d1]->p_lk_left;
1953	}
1954	}
1955
1956	/*
1957	** Temporarily set the t_edge lengths and update transition prob's at the
1958	** prune site.
1959	*/
1960	v_prune->b[d1]->l_old->v = v_prune->b[d1]->l->v;
1961	v_prune->b[d2]->l_old->v = v_prune->b[d2]->l->v;
1962	v_prune->b[d1]->l->v = l_connect;
1963	v_prune->b[d2]->l->v = l_connect;
1964	Update_PMat_At_Given_Edge (v_prune->b[d1], tree);
1965	Update_PMat_At_Given_Edge (v_prune->b[d2], tree);
1966
1967	/*
1968	** Get the relevant average subtree distance within the pruned subtree.
1969	*/
1970	if (!u_prune->tax)
1971	{
1972	u1 = u2 = NULL;
1973	for (i = 0; i < 3; i++)
1974	{
1975	if (u_prune->b[i] != e_prune)
1976	{
1977	if (u1 == NULL)
1978	{
1979	u1 = u_prune->v[i];
1980	}
1981	else
1982	{
1983	u2 = u_prune->v[i];
1984	}
1985	}
1986	}
1987	d_uu = subtree_dist[u1->num][u2->num];
1988	}
1989	else
1990	{
1991	d_uu = 0.0;
1992	}
1993
1994	/*
1995	** Try the best candidate SPRs and estimate the change in likelihood.
1996	*/
1997	best_d_lk = -1.0*BIG;
1998	best_cand = 0;
1999	for (cand = 0; cand < tree->mod->s_opt->wim_n_best; cand++)
2000	{
2001	/*
2002	** If there are no more candidates, bail out...
2003	*/
2004	if (FABS(rgrft_cand[cand]->d_L - 1.0*BIG) < SMALL)
2005	{
2006	break;
2007	}
2008	else
2009	{
2010	nr_d_lk++;
2011	}
2012
2013	/*
2014	** Get the relevant nodes and edges.
2015	*/
2016	v_n = rgrft_cand[cand]->v_n;
2017	v_nx1 = rgrft_cand[cand]->v_nx1;
2018	e_regraft = rgrft_cand[cand]->e_regraft;
2019
2020	/*
2021	** Update the relevant partial likelihoods along the path between the prune
2022	** and regraft positions (temporarily save the first one).
2023	*/
2024	n = rgrft_cand[cand]->dist;
2025	e_tmp = NULL;
2026	p_lk = NULL;
2027	p_sum = NULL;
2028	for (i = 1; i <= n; i++)
2029	{
2030	/*
2031	** Get the next t_edge along the path.
2032	*/
2033	for (j = 0; j < 3; j++)
2034	{
2035	if (rgrft_cand[cand]->path[i]->v[j] == rgrft_cand[cand]->path[i+1])
2036	{
2037	e_tmp = rgrft_cand[cand]->path[i]->b[j];
2038	break;
2039	}
2040	}
2041	if (i == 1)
2042	{
2043	/*
2044	** Save the first partial likelihood along the path.
2045	*/
2046	if (rgrft_cand[cand]->path[i] == e_tmp->left)
2047	{
2048	p_lk = e_tmp->p_lk_left;
2049	p_sum = e_tmp->sum_scale_left;
2050	}
2051	else
2052	{
2053	p_lk = e_tmp->p_lk_rght;
2054	p_sum = e_tmp->sum_scale_rght;
2055	}
2056	for (pat = 0; pat < tree->n_pattern; pat++)
2057	{
2058	sum_scale_tmp[pat] = p_sum[pat];
2059	for (cat = 0; cat < tree->mod->ras->n_catg; cat++)
2060	{
2061	for (ste = 0; ste < tree->mod->ns; ste++)
2062	{
2063	p_lk_tmp[patdim1+catdim2+ste] = p_lk[patdim1+catdim2+ste];
2064	}
2065	}
2066	}
2067	}
2068	Update_P_Lk (tree, e_tmp, rgrft_cand[cand]->path[i]);
2069	}
2070	if (v_n == e_regraft->left)
2071	{
2072	v_tmp->b[1]->p_lk_rght = e_regraft->p_lk_left;
2073	v_tmp->b[2]->p_lk_rght = e_regraft->p_lk_rght;
2074	v_tmp->b[1]->sum_scale_rght = e_regraft->sum_scale_left;
2075	v_tmp->b[2]->sum_scale_rght = e_regraft->sum_scale_rght;
2076	}
2077	else
2078	{
2079	v_tmp->b[1]->p_lk_rght = e_regraft->p_lk_rght;
2080	v_tmp->b[2]->p_lk_rght = e_regraft->p_lk_left;
2081	v_tmp->b[1]->sum_scale_rght = e_regraft->sum_scale_rght;
2082	v_tmp->b[2]->sum_scale_rght = e_regraft->sum_scale_left;
2083	}
2084
2085	/*
2086	** Estimate t_edge lengths of the three relevant regraft edges based on
2087	** average subtree distances.
2088	**
2089	** l_01
2090	*/
2091	/*
2092	** Alternative method of estimating l_01. Kept it around for reference...
2093	**
2094	l_01 = subtree_dist[u_prune->num][u_n->num] - (0.5 * d_uu);
2095	if (!u_n->tax)
2096	{
2097	u1 = u2 = NULL;
2098	for (i = 0; i < 3; i++)
2099	{
2100	if (u_n->v[i] != v_n)
2101	{
2102	if (u1 == NULL)
2103	{
2104	u1 = u_n->v[i];
2105	}
2106	else
2107	{
2108	u2 = u_n->v[i];
2109	}
2110	}
2111	}
2112	l_01 -= 0.5 * subtree_dist[u1->num][u2->num];
2113	}
2114	for (i = 0; i < 3; i++)
2115	{
2116	if (u_n->v[i] == v_n)
2117	{
2118	l_01 -= u_n->b[i]->l->v;
2119	break;
2120	}
2121	}
2122	*/
2123	l_01 = rgrft_cand[cand]->d_up_v - (0.5 * rgrft_cand[cand]->d_un_v) -
2124	(0.5 * d_uu);
2125	/*
2126	** l_02
2127	*/
2128	l_02 = subtree_dist[u_prune->num][v_nx1->num] - (0.5 * d_uu);
2129	if (!v_nx1->tax)
2130	{
2131	u1 = u2 = NULL;
2132	for (i = 0; i < 3; i++)
2133	{
2134	if (v_nx1->v[i] != v_n)
2135	{
2136	if (u1 == NULL)
2137	{
2138	u1 = v_nx1->v[i];
2139	}
2140	else
2141	{
2142	u2 = v_nx1->v[i];
2143	}
2144	}
2145	}
2146	l_02 -= (0.5 * subtree_dist[u1->num][u2->num]);
2147	}
2148	/*
2149	** l_12
2150	*/
2151	l_12 = e_regraft->l->v;
2152	/*
2153	** Simple estimates.
2154	*/
2155	l_simple[0] = l_02 - (0.5*e_regraft->l->v);
2156	l_simple[1] = 0.5 * e_regraft->l->v;
2157	l_simple[2] = 0.5 * e_regraft->l->v;
2158	/*
2159	** Average subtree distance based estimates.
2160	*/
2161	l_dist[0] = 0.5 * ( l_01 + l_02 - l_12);
2162	l_dist[1] = 0.5 * ( l_01 - l_02 + l_12);
2163	l_dist[2] = 0.5 * (-l_01 + l_02 + l_12);
2164	/*
2165	** Take the average of the two estimates.
2166	*/
2167	l_est[0] = (l_simple[0] + l_dist[0]) / 2.0;
2168	l_est[1] = (l_simple[1] + l_dist[1]) / 2.0;
2169	l_est[2] = (l_simple[2] + l_dist[2]) / 2.0;
2170
2171	/*
2172	** Set the t_edge lengths and update the relevant transition prob's and
2173	** partial likelihoods in the temporary regraft structure.
2174	*/
2175	for (i = 0; i < 3; i++)
2176	{
2177	v_tmp->b[i]->l->v = l_est[i];
2178	Update_PMat_At_Given_Edge (v_tmp->b[i], tree);
2179	}
2180	Update_P_Lk (tree, v_tmp->b[0], v_tmp);
2181
2182	/*
2183	** Calculate the change in likelihood locally. Save it and the estimated edge
2184	** lengths in the current candidate in the list.
2185	*/
2186	new_lk = Lk(v_tmp->b[0],tree);
2187	rgrft_cand[cand]->delta_lk = new_lk - cur_lk;
2188	rgrft_cand[cand]->rgrft_rank = cand;
2189	rgrft_cand[cand]->optim_rank = -1;
2190	rgrft_cand[cand]->globl_rank = -1;
2191	rgrft_cand[cand]->l_connect = l_connect;
2192	for (i = 0; i < 3; i++)
2193	{
2194	rgrft_cand[cand]->l_est[i] = v_tmp->b[i]->l->v;
2195	}
2196	if (rgrft_cand[cand]->delta_lk > best_d_lk)
2197	{
2198	best_d_lk = rgrft_cand[cand]->delta_lk;
2199	best_cand = cand;
2200	}
2201
2202	/*
2203	** Reset the partial likelihoods along the path.
2204	*/
2205	for (pat = 0; pat < tree->n_pattern; pat++)
2206	{
2207	p_sum[pat] = sum_scale_tmp[pat];
2208	for (cat = 0; cat < tree->mod->ras->n_catg; cat++)
2209	{
2210	for (ste = 0; ste < tree->mod->ns; ste++)
2211	{
2212	p_lk[patdim1+catdim2+ste] = p_lk_tmp[patdim1+catdim2+ste];
2213	}
2214	}
2215	}
2216	n = rgrft_cand[cand]->dist;
2217	for (i = 2; i <= n; i++)
2218	{
2219	for (j = 0; j < 3; j++)
2220	{
2221	if (rgrft_cand[cand]->path[i]->v[j] == rgrft_cand[cand]->path[i+1])
2222	{
2223	e_tmp = rgrft_cand[cand]->path[i]->b[j];
2224	break;
2225	}
2226	}
2227	Update_P_Lk (tree, e_tmp, rgrft_cand[cand]->path[i]);
2228	}
2229	}
2230
2231	/*
2232	** If the best candidate is within the tree->mod->s_opt->wim_n_optim best ones, save it in the list of
2233	** optimization candidates.
2234	*/
2235	if (rgrft_cand[best_cand]->delta_lk > optim_cand[tree->mod->s_opt->wim_n_optim-1]->delta_lk)
2236	{
2237	i = tree->mod->s_opt->wim_n_optim-1;
2238	optim_cand[i]->v_prune = rgrft_cand[best_cand]->v_prune;
2239	optim_cand[i]->u_prune = rgrft_cand[best_cand]->u_prune;
2240	optim_cand[i]->v_n = rgrft_cand[best_cand]->v_n;
2241	optim_cand[i]->v_nx1 = rgrft_cand[best_cand]->v_nx1;
2242	optim_cand[i]->u_n = rgrft_cand[best_cand]->u_n;
2243	optim_cand[i]->e_prune = rgrft_cand[best_cand]->e_prune;
2244	optim_cand[i]->e_regraft = rgrft_cand[best_cand]->e_regraft;
2245	optim_cand[i]->d_L = rgrft_cand[best_cand]->d_L;
2246	optim_cand[i]->dist = rgrft_cand[best_cand]->dist;
2247	optim_cand[i]->rgrft_rank = rgrft_cand[best_cand]->rgrft_rank;
2248	optim_cand[i]->optim_rank = rgrft_cand[best_cand]->optim_rank;
2249	optim_cand[i]->globl_rank = rgrft_cand[best_cand]->globl_rank;
2250	optim_cand[i]->l_connect = rgrft_cand[best_cand]->l_connect;
2251
2252	for (j = 0; j < 3; j++)
2253	{
2254	optim_cand[i]->l_est[j] = rgrft_cand[best_cand]->l_est[j];
2255	}
2256	optim_cand[i]->delta_lk = rgrft_cand[best_cand]->delta_lk;
2257	/*
2258	** Move the candidate to the appropriate position in the list, so the list
2259	** remains sorted in decreasing delta_Lk value.
2260	*/
2261	while ((i > 0) && (optim_cand[i]->delta_lk > optim_cand[i-1]->delta_lk))
2262	{
2263	tmp_cand = optim_cand[i];
2264	optim_cand[i] = optim_cand[i-1];
2265	optim_cand[i-1] = tmp_cand;
2266	i--;
2267	}
2268	}
2269
2270	/*
2271	** Swap back the relevant partial likelihoods at the prune site.
2272	*/
2273	if (v_prune == v_prune->b[d1]->left)
2274	{
2275	v_prune->b[d1]->p_lk_left = p_lk1_tmp;
2276	}
2277	else
2278	{
2279	v_prune->b[d1]->p_lk_rght = p_lk1_tmp;
2280	}
2281	if (v_prune == v_prune->b[d2]->left)
2282	{
2283	v_prune->b[d2]->p_lk_left = p_lk2_tmp;
2284	}
2285	else
2286	{
2287	v_prune->b[d2]->p_lk_rght = p_lk2_tmp;
2288	}
2289
2290	/*
2291	** Reset the relevant t_edge lengths and transition prob's at the prune site.
2292	*/
2293	v_prune->b[d1]->l->v = v_prune->b[d1]->l_old->v;
2294	v_prune->b[d2]->l->v = v_prune->b[d2]->l_old->v;
2295	Update_PMat_At_Given_Edge (v_prune->b[d1], tree);
2296	Update_PMat_At_Given_Edge (v_prune->b[d2], tree);
2297
2298	/*
2299	** Return the best candidate.
2300	*/
2301	return (best_cand);
2302	}
2303
2304
2305	/*
2306	** Make_Move: Perform an actual SPR move and calculate the new likelihood.
2307	**
2308	** Parameters:
2309	** - candidate: The candidate move to perform.
2310	** - tree: The tree on which to perform the move.
2311	**
2312	*/
2313
2314	void Make_Move (_move_ move, int type, t_tree tree)
2315	{
2316	int i;
2317	t_node v_prune, u_prune, v_n, root;
2318	t_edge e_prune, e_regraft, e_connect, e_avail;
2319	phydbl new_lk;
2320
2321	/*
2322	** Get the relevant nodes and edges.
2323	*/
2324	v_prune = move->v_prune;
2325	u_prune = move->u_prune;
2326	v_n = move->v_n;
2327	e_prune = move->e_prune;
2328	e_regraft = move->e_regraft;
2329	/* PhyML_Printf (" making move: %d -> %d (%f)\n", e_prune->num, e_regraft->num, move->delta_lk); */
2330	/*
2331	** Perform the move and set t_edge lengths.
2332	*/
2333	Prune (e_prune, v_prune, &(e_connect), &(e_avail), tree);
2334	Regraft (e_regraft, v_prune, e_avail, tree);
2335	e_connect->l->v = move->l_connect;
2336
2337	for (i = 0; i < 3; i++)
2338	{
2339	if (v_prune->v[i] == u_prune)
2340	{
2341	v_prune->b[i]->l->v = move->l_est[0];
2342	}
2343	else if (v_prune->v[i] == v_n)
2344	{
2345	v_prune->b[i]->l->v = move->l_est[1];
2346	}
2347	else
2348	{
2349	v_prune->b[i]->l->v = move->l_est[2];
2350	}
2351	}
2352
2353
2354	if(type > 0) /* local or global move */
2355	{
2356	Restore_Br_Len(tree);
2357	}
2358
2359	/*
2360	** Calculate the new likelihood.
2361	*/
2362	Set_Both_Sides(YES,tree);
2363	new_lk = Lk(NULL,tree);
2364
2365	if(tree->c_lnL < cur_lk-tree->mod->s_opt->min_diff_lk_local)
2366	{
2367	PhyML_Printf("\n. tree->c_lnL = %f cur_lk = %f",tree->c_lnL,cur_lk);
2368	PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
2369	Warn_And_Exit("");
2370	}
2371	cur_lk = new_lk;
2372
2373	/*
2374	** Recalculate the average distances between all (non-overlapping) subtrees.
2375	*/
2376	root = tree->a_nodes[0];
2377	PostOrder_v (tree, root->v[2], root->b[2]);
2378	}
2379
2380
2381	/*
2382	** Find_Optim_Local: Perform local t_edge length optimization on the candidates in the
2383	** optimization list, and return the first one that gives an
2384	** improvement in likelihood.
2385	**
2386	** Parameters:
2387	** - tree: The tree on which to check the moves.
2388	**
2389	** Returns:
2390	** If an improvement was found: The candidate that gives the improvement.
2391	** Otherwise: -1.
2392	*/
2393
2394	int Find_Optim_Local (t_tree *tree)
2395	{
2396	int best_cand, cand, i;
2397	t_node v_prune, u_prune, *v_n;
2398	t_edge e_prune, e_regraft, e_connect, e_avail;
2399	phydbl max_change, new_lk;
2400	_move_ move, tmp_cand;
2401
2402	/*
2403	** Try all candidate moves starting from the first one.
2404	*/
2405	best_cand = -1;
2406	max_change = 1.0/BIG;
2407	for(cand = 0; cand < tree->mod->s_opt->wim_n_optim; cand++)
2408	{
2409	move = optim_cand[cand];
2410	if(move->delta_lk > -1.0*BIG)
2411	{
2412	/*
2413	** Get the relevant nodes and edges.
2414	*/
2415	v_prune = move->v_prune;
2416	u_prune = move->u_prune;
2417	v_n = move->v_n;
2418	e_prune = move->e_prune;
2419	e_regraft = move->e_regraft;
2420
2421	/*
2422	** Perform the move and set t_edge lengths.
2423	*/
2424	Prune (e_prune, v_prune, &(e_connect), &(e_avail), tree);
2425	Regraft (e_regraft, v_prune, e_avail, tree);
2426	e_connect->l_old->v = e_connect->l->v;
2427	e_connect->l->v = move->l_connect;
2428
2429	for (i = 0; i < 3; i++)
2430	{
2431	v_prune->b[i]->l_old->v = v_prune->b[i]->l->v;
2432
2433	if (v_prune->v[i] == u_prune)
2434	{
2435	v_prune->b[i]->l->v = move->l_est[0];
2436	}
2437	else if (v_prune->v[i] == v_n)
2438	{
2439	v_prune->b[i]->l->v = move->l_est[1];
2440	}
2441	else
2442	{
2443	v_prune->b[i]->l->v = move->l_est[2];
2444	}
2445	}
2446
2447	Set_Both_Sides(YES,tree);
2448	Lk(NULL,tree); // Not sure anymore whether this is required...
2449
2450	/*
2451	** Use Brent optimization on the relevant edges at the regraft position
2452	** and calculate the new likelihood value.
2453	*/
2454	Br_Len_Brent (0.05,20.,v_prune->b[0], tree);
2455	Br_Len_Brent (0.05,20.,v_prune->b[1], tree);
2456	Br_Len_Brent (0.05,20.,v_prune->b[2], tree);
2457
2458	/* Update_PMat_At_Given_Edge (v_prune->b[0], tree); */
2459	/* Update_PMat_At_Given_Edge (v_prune->b[1], tree); */
2460	/* Update_PMat_At_Given_Edge (v_prune->b[2], tree); */
2461
2462	Update_P_Lk (tree, v_prune->b[0], v_prune);
2463	new_lk = Lk(v_prune->b[0],tree);
2464
2465	/* PhyML_Printf("\n. local new_lk = %f",new_lk); */
2466
2467	/*
2468	** Save the change in likelihood and move the current candidate to the
2469	** appropriate place in the list.
2470	*/
2471	move->delta_lk = new_lk - cur_lk;
2472	move->optim_rank = cand;
2473	i = cand;
2474	while ((i > 0) && (optim_cand[i]->delta_lk > optim_cand[i-1]->delta_lk))
2475	{
2476	tmp_cand = optim_cand[i];
2477	optim_cand[i] = optim_cand[i-1];
2478	optim_cand[i-1] = tmp_cand;
2479	i--;
2480	}
2481	if (move->delta_lk > max_change)
2482	{
2483	best_cand = i;
2484	max_change = move->delta_lk;
2485	Record_Br_Len(tree);
2486	}
2487
2488
2489	/*
2490	** Undo the move again.
2491	*/
2492	Prune (e_prune, v_prune, &(e_regraft), &(e_avail), tree);
2493	Regraft (e_connect, v_prune, e_avail, tree);
2494	e_regraft->l->v = e_regraft->l_old->v;
2495	for (i = 0; i < 3; i++)
2496	{
2497	v_prune->b[i]->l->v = v_prune->b[i]->l_old->v;
2498	}
2499	Set_Both_Sides(YES,tree);
2500	Lk(NULL,tree);
2501	nr_loc++;
2502	/* PhyML_Printf("\n. local back to = %f",tree->c_lnL); */
2503	}
2504	else
2505	{
2506	break;
2507	}
2508
2509	/*
2510	** If an improvement was found, forget the other candidates...
2511	*/
2512	if (best_cand >= 0)
2513	{
2514	break;
2515	}
2516	}
2517
2518	/*
2519	** Return the best candidate.
2520	*/
2521	return (best_cand);
2522	}
2523
2524
2525	/*
2526	** Find_Optim_Globl: Perform global t_edge length optimization on the candidates in the
2527	** optimization list, and return the first one that gives an
2528	** improvement in likelihood.
2529	**
2530	** Parameters:
2531	** - tree: The tree on which to check the moves.
2532	**
2533	** Returns:
2534	** If an improvement is found: The candidate that gives the improvement.
2535	** Otherwise: -1.
2536	*/
2537
2538	int Find_Optim_Globl (t_tree *tree)
2539	{
2540	int best_cand, cand, i;
2541	t_node v_prune, u_prune, *v_n;
2542	t_edge e_prune, e_regraft, e_connect, e_avail;
2543	phydbl max_change, new_lk;
2544	_move_ *move;
2545
2546	/*
2547	** Try all moves.
2548	*/
2549	best_cand = -1;
2550	max_change = 1.0/BIG;
2551	for (cand = 0; cand < tree->mod->s_opt->wim_n_globl; cand++)
2552	{
2553	move = optim_cand[cand];
2554	if (move->delta_lk > -1.0*BIG)
2555	{
2556	Record_Br_Len(tree);
2557
2558	/*
2559	** Get the relevant nodes and edges.
2560	*/
2561	v_prune = move->v_prune;
2562	u_prune = move->u_prune;
2563	v_n = move->v_n;
2564	e_prune = move->e_prune;
2565	e_regraft = move->e_regraft;
2566
2567	/*
2568	** Perform the move and optimize all t_edge lengths.
2569	*/
2570	Prune (e_prune, v_prune, &(e_connect), &(e_avail), tree);
2571	Regraft (e_regraft, v_prune, e_avail, tree);
2572	e_connect->l_old->v = e_connect->l->v;
2573	e_connect->l->v = move->l_connect;
2574
2575	for (i = 0; i < 3; i++)
2576	{
2577	v_prune->b[i]->l_old->v = v_prune->b[i]->l->v;
2578	if (v_prune->v[i] == u_prune)
2579	{
2580	v_prune->b[i]->l->v = move->l_est[0];
2581	}
2582	else if (v_prune->v[i] == v_n)
2583	{
2584	v_prune->b[i]->l->v = move->l_est[1];
2585	}
2586	else
2587	{
2588	v_prune->b[i]->l->v = move->l_est[2];
2589	}
2590	}
2591
2592	Set_Both_Sides(YES,tree);
2593	Lk(NULL,tree);
2594	Optimize_Br_Len_Serie (tree);
2595	Set_Both_Sides(YES,tree);
2596	new_lk = Lk (NULL,tree);
2597
2598	/* PhyML_Printf("\n. global new_lk = %f\n",tree->c_lnL); */
2599
2600	/*
2601	** Save the change in likelihood and undo the move.
2602	*/
2603
2604	move->delta_lk = new_lk - cur_lk;
2605	move->globl_rank = cand;
2606	if (move->delta_lk > max_change)
2607	{
2608	best_cand = cand;
2609	max_change = move->delta_lk;
2610	Record_Br_Len(tree);
2611	}
2612
2613	Prune (e_prune, v_prune, &(e_regraft), &(e_avail), tree);
2614	Regraft (e_connect, v_prune, e_avail, tree);
2615	e_regraft->l->v = e_regraft->l_old->v;
2616	for (i = 0; i < 3; i++)
2617	{
2618	v_prune->b[i]->l->v = v_prune->b[i]->l_old->v;
2619	}
2620	Set_Both_Sides(YES,tree);
2621	Restore_Br_Len(tree);
2622	Lk(NULL,tree);
2623	nr_glb++;
2624	/* PhyML_Printf("\n. global back to = %f",tree->c_lnL); */
2625
2626	}
2627	else break;
2628
2629	/*
2630	** If an improvement was found, forget the other candidates...
2631	*/
2632	if (best_cand >= 0) break;
2633	}
2634
2635	/*
2636	** Return the best candidate.
2637	*/
2638	return (best_cand);
2639	}
2640
2641
2642	/*
2643	** Prune: Prune the subtree at a certain t_edge and node. Note that edge
2644	** lengths are not set and partial likelihoods are not updated.
2645	**
2646	** Parameters:
2647	** - e: The t_edge at which to prune the subtree.
2648	** - v: The t_node adjacent to t_edge e which forms the root of the subtree.
2649	** - e_connect: An t_edge pointer which will point to the t_edge that was left
2650	** after pruning.
2651	** - e_avail: The t_edge that is "left over" and should be used in the
2652	** regrafting step.
2653	** - tree: The tree on which the pruning is done.
2654	**
2655	**
2656	**
2657	** \ /
2658	** o
2659	** \|
2660	** \| e --> subtree to be pruned
2661	** \|
2662	** o v
2663	** / \
2664	** e1/ \e2
2665	** / \
2666	** o o
2667	** u1 u2 --> such that u1->num < u2->num
2668	*/
2669
2670	void Prune (t_edge e, t_node v, t_edge e_connect, t_edge e_avail, t_tree *tree)
2671	{
2672	int dir0, dir1, dir2, v0, v1, v2, tmp_dir, i, j, k;
2673	t_node u1, u2, *tmp_node;
2674	t_edge e1, e2;
2675	int *sum_scale_f;
2676	phydbl *p_lk;
2677	int dim1, dim2;
2678
2679
2680	dim1 = tree->mod->ns * tree->mod->ras->n_catg;
2681	dim2 = tree->mod->ns;
2682
2683	/*
2684	** Get the relevant directions, nodes and edges.
2685	** Make sure that t_node u1 is the t_node with the smaller number.
2686	*/
2687	dir0 = -1;
2688	for (i = 0; i < 3; i++)
2689	{
2690	if (v->b[i] == e)
2691	{
2692	dir0 = i;
2693	break;
2694	}
2695	}
2696	dir1 = (dir0 + 1) % 3;
2697	dir2 = 3 - dir0 - dir1;
2698	u1 = v->v[dir1];
2699	u2 = v->v[dir2];
2700	if (u1->num > u2->num)
2701	{
2702	tmp_node = u1;
2703	u1 = u2;
2704	u2 = tmp_node;
2705	tmp_dir = dir1;
2706	dir1 = dir2;
2707	dir2 = tmp_dir;
2708	}
2709	e1 = v->b[dir1];
2710	e2 = v->b[dir2];
2711
2712	/*
2713	** Detach t_node v from the tree.
2714	*/
2715	v->v[dir1] = NULL;
2716	v->v[dir2] = NULL;
2717	v->b[dir1] = NULL;
2718	v->b[dir2] = NULL;
2719
2720	/*
2721	** Connect nodes u1 and u2 via t_edge e1 and copy relevant partial likelihoods.
2722	*/
2723	if (u2 == e2->left)
2724	{
2725	v0 = e2->l_r;
2726	v1 = e2->l_v1;
2727	v2 = e2->l_v2;
2728	sum_scale_f = e2->sum_scale_left;
2729	p_lk = e2->p_lk_left;
2730	}
2731	else
2732	{
2733	v0 = e2->r_l;
2734	v1 = e2->r_v1;
2735	v2 = e2->r_v2;
2736	sum_scale_f = e2->sum_scale_rght;
2737	p_lk = e2->p_lk_rght;
2738	}
2739	if (u1 == e1->left)
2740	{
2741	e1->rght = u2;
2742	e1->r_l = v0;
2743	e1->r_v1 = v1;
2744	e1->r_v2 = v2;
2745	for (i = 0; i < tree->n_pattern; i++)
2746	{
2747	e1->sum_scale_rght[i] = sum_scale_f[i];
2748	for (j = 0; j < tree->mod->ras->n_catg; j++)
2749	{
2750	for (k = 0; k < tree->mod->ns; k++)
2751	{
2752	e1->p_lk_rght[idim1+jdim2+k] = p_lk[idim1+jdim2+k];
2753	}
2754	}
2755	}
2756	}
2757	else
2758	{
2759	e1->left = u2;
2760	e1->l_r = v0;
2761	e1->l_v1 = v1;
2762	e1->l_v2 = v2;
2763	for (i = 0; i < tree->n_pattern; i++)
2764	{
2765	e1->sum_scale_left[i] = sum_scale_f[i];
2766	for (j = 0; j < tree->mod->ras->n_catg; j++)
2767	{
2768	for (k = 0; k < tree->mod->ns; k++)
2769	{
2770	e1->p_lk_left[idim1+jdim2+k] = p_lk[idim1+jdim2+k];
2771	}
2772	}
2773	}
2774	}
2775	for (i = 0; i < 3; i++)
2776	{
2777	if (u1->v[i] == v)
2778	{
2779	u1->v[i] = u2;
2780	}
2781	if (u2->v[i] == v)
2782	{
2783	u2->v[i] = u1;
2784	u2->b[i] = e1;
2785	u2->l[i] = e1->l->v;
2786	}
2787	}
2788
2789	/*
2790	** Make sure that a possible tip t_node is still on the right side.
2791	*/
2792	if (e1->left->tax)
2793	{
2794	/*
2795	** Swap left and right.
2796	*/
2797	tmp_node = e1->left;
2798	e1->left = e1->rght;
2799	e1->rght = tmp_node;
2800	tmp_dir = e1->l_r;
2801	e1->l_r = e1->r_l;
2802	e1->r_l = tmp_dir;
2803	tmp_dir = e1->l_v1;
2804	e1->l_v1 = e1->r_v1;
2805	e1->r_v1 = tmp_dir;
2806	tmp_dir = e1->l_v2;
2807	e1->l_v2 = e1->r_v2;
2808	e1->r_v2 = tmp_dir;
2809	p_lk = e1->p_lk_left;
2810	e1->p_lk_left = e1->p_lk_rght;
2811	e1->p_lk_rght = p_lk;
2812	sum_scale_f = e1->sum_scale_left;
2813	e1->sum_scale_left = e1->sum_scale_rght;
2814	e1->sum_scale_rght = sum_scale_f;
2815	}
2816
2817	/*
2818	** Set the connecting and available edges.
2819	*/
2820	*(e_connect) = e1;
2821	*(e_avail) = e2;
2822	}
2823
2824
2825	/*
2826	** Regraft: Regraft a subtree at a certain edge. Note that t_edge lengths
2827	** are not set and partial likelihoods are not updated.
2828	**
2829	** Parameters:
2830	** - e: The t_edge to regraft the subtree on.
2831	** - v: The root of the subtree to regraft.
2832	** - avail: A previously deleted t_edge now available for insertion again.
2833	** - tree: The tree on which the regrafting is done.
2834	**
2835	**
2836	** \ /
2837	** o
2838	** \|
2839	** \| --> subtree to regraft
2840	** \|
2841	** o v
2842	**
2843	** o--------o
2844	** u1 e u2 --> such that u1->num < u2->num
2845	*/
2846
2847	void Regraft (t_edge e, t_node v, t_edge avail, t_tree tree)
2848	{
2849	int dir0, dir1, dir2, i, j, k;
2850	int *sum_scale_f;
2851	phydbl *p_lk;
2852	t_node u1, u2;
2853	int dim1, dim2;
2854
2855	dim1 = tree->mod->ns * tree->mod->ras->n_catg;
2856	dim2 = tree->mod->ns ;
2857
2858	/*
2859	** Get the relevant directions and nodes.
2860	*/
2861	dir0 = -1;
2862	for (i = 0; i < 3; i++)
2863	{
2864	if (v->b[i] != NULL)
2865	{
2866	dir0 = i;
2867	break;
2868	}
2869	}
2870	dir1 = (dir0 + 1) % 3;
2871	dir2 = 3 - dir0 - dir1;
2872	if (e->left->num < e->rght->num)
2873	{
2874	u1 = e->left;
2875	u2 = e->rght;
2876	sum_scale_f = e->sum_scale_rght;
2877	p_lk = e->p_lk_rght;
2878	}
2879	else
2880	{
2881	u1 = e->rght;
2882	u2 = e->left;
2883	sum_scale_f = e->sum_scale_left;
2884	p_lk = e->p_lk_left;
2885	}
2886
2887	/*
2888	** Connect nodes v and u2 via the available t_edge 'avail' and copy the
2889	** relevant partial likelihood.
2890	** (We want to do this first, cos we need some of the values of edge
2891	** e before changing them below).
2892	*/
2893	avail->left = v;
2894	avail->rght = u2;
2895	avail->l_r = dir2;
2896	avail->l_v1 = dir0;
2897	avail->l_v2 = dir1;
2898	v->v[dir2] = u2;
2899	v->b[dir2] = avail;
2900	for (i = 0; i < 3; i++)
2901	{
2902	if (e == u2->b[i])
2903	{
2904	u2->v[i] = v;
2905	u2->b[i] = avail;
2906	avail->r_l = i;
2907	avail->r_v1 = (i + 1) % 3;
2908	avail->r_v2 = 3 - i - avail->r_v1;
2909	break;
2910	}
2911	}
2912	for (i = 0; i < tree->n_pattern; i++)
2913	{
2914	avail->sum_scale_rght[i] = sum_scale_f[i];
2915	for (j = 0; j < tree->mod->ras->n_catg; j++)
2916	{
2917	for (k = 0; k < tree->mod->ns; k++)
2918	{
2919	avail->p_lk_rght[idim1+jdim2+k] = p_lk[idim1+jdim2+k];
2920	}
2921	}
2922	}
2923
2924	/*
2925	** Connect nodes v and u1 via t_edge e.
2926	*/
2927	if (u1 == e->left)
2928	{
2929	e->rght = v;
2930	e->r_l = dir1;
2931	e->r_v1 = dir0;
2932	e->r_v2 = dir2;
2933	}
2934	else
2935	{
2936	e->left = v;
2937	e->l_r = dir1;
2938	e->l_v1 = dir0;
2939	e->l_v2 = dir2;
2940	}
2941	v->v[dir1] = u1;
2942	v->b[dir1] = e;
2943	for (i = 0; i < 3; i++)
2944	{
2945	if (u1->v[i] == u2)
2946	{
2947	u1->v[i] = v;
2948	break;
2949	}
2950	}
2951	}
2952
2953
2954	/*
2955	** PostOrder_v: Recursively visit all nodes v in postorder to calculate
2956	** the average distance between subtrees.
2957	**
2958	** Parameters:
2959	** - tree: The tree for which to calculate the average distances.
2960	** - v: The current node.
2961	** - e: The t_edge we came from.
2962	*/
2963
2964	void PostOrder_v (t_tree tree, t_node v, t_edge *e)
2965	{
2966	int i;
2967	t_node *w;
2968
2969	/*
2970	** If v is not a taxon, recurse.
2971	*/
2972	if (!v->tax)
2973	{
2974	for (i = 0; i < 3; i++)
2975	{
2976	if (v->b[i] != e)
2977	{
2978	PostOrder_v (tree, v->v[i], v->b[i]);
2979	}
2980	}
2981	}
2982
2983	/*
2984	** Recurse over all nodes w not in the current subtree and calculate
2985	** the average distance between the current subtree and all others.
2986	*/
2987	if (v == e->left)
2988	{
2989	w = e->rght;
2990	}
2991	else
2992	{
2993	w = e->left;
2994	}
2995	PostOrder_w (tree, v, e, w, e);
2996	}
2997
2998
2999	/*
3000	** PostOrder_w: Recursively visit all nodes w not in the subtree of v in
3001	** postorder and calculate the average distance between the
3002	** subtree of v and all others.
3003	**
3004	** Parameters:
3005	** - tree: The tree for which to calculate the average distances.
3006	** - v: The root t_node of the first subtree.
3007	** - v_e: The t_edge adjacent to the first subtree.
3008	** - w: The current node.
3009	** - e: The t_edge we came from.
3010	*/
3011
3012	void PostOrder_w (t_tree tree, t_node v, t_edge v_e, t_node w, t_edge *e)
3013	{
3014	int i;
3015	t_node w1, w2, v1, v2;
3016
3017	/*
3018	** If w is not a taxon, recurse.
3019	*/
3020	if (!w->tax)
3021	{
3022	for (i = 0; i < 3; i++)
3023	{
3024	if (w->b[i] != e)
3025	{
3026	PostOrder_w (tree, v, v_e, w->v[i], w->b[i]);
3027	}
3028	}
3029	}
3030
3031	/*
3032	** Calculate the average distance between the subtrees defined by
3033	** nodes v and w.
3034	*/
3035	if (v->tax && w->tax)
3036	{
3037	subtree_dist[v->num][w->num] = seq_dist->dist[v->num][w->num];
3038	}
3039	else if (!v->tax)
3040	{
3041	v1 = v2 = NULL;
3042	for (i = 0; i < 3; i++)
3043	{
3044	if (v->b[i] != v_e)
3045	{
3046	if (v1 == NULL)
3047	{
3048	v1 = v->v[i];
3049	}
3050	else
3051	{
3052	v2 = v->v[i];
3053	}
3054	}
3055	}
3056	subtree_dist[v->num][w->num] = 0.5*(subtree_dist[v1->num][w->num] +
3057	subtree_dist[v2->num][w->num]);
3058	}
3059	else
3060	{
3061	w1 = w2 = NULL;
3062	for (i = 0; i < 3; i++)
3063	{
3064	if (w->b[i] != e)
3065	{
3066	if (w1 == NULL)
3067	{
3068	w1 = w->v[i];
3069	}
3070	else
3071	{
3072	w2 = w->v[i];
3073	}
3074	}
3075	}
3076	subtree_dist[v->num][w->num] = 0.5*(subtree_dist[v->num][w1->num] +
3077	subtree_dist[v->num][w2->num]);
3078	}
3079	subtree_dist[w->num][v->num] = subtree_dist[v->num][w->num];
3080	}
3081
3082
3083	/*********************************************************/
3084	/*********************************************************/
3085	/*********************************************************/
3086	/* Below are my functions for SPR search (Stephane Guindon, 2007) */
3087
3088
3089	void Randomize_Spr_List(t_tree *tree)
3090	{
3091	int i,j;
3092	t_spr *buff;
3093
3094	For(i,tree->size_spr_list)
3095	{
3096	j = (int)FLOOR(rand()/(RAND_MAX+1.)*tree->size_spr_list);
3097	buff = tree->spr_list[i];
3098	tree->spr_list[i] = tree->spr_list[j];
3099	tree->spr_list[j] = buff;
3100	}
3101	}
3102
3103	/*********************************************************/
3104
3105	int Spr(phydbl init_lnL, t_tree *tree)
3106	{
3107	int br;
3108	int pars_diff, max_pars_diff, new_pars, old_pars;
3109	t_edge *b;
3110
3111	Set_Both_Sides(YES,tree);
3112	pars_diff = -1;
3113	max_pars_diff = -1;
3114	new_pars = -1;
3115	old_pars = -1;
3116
3117	Reset_Spr_List(tree);
3118
3119	For(br,2*tree->n_otu-3)
3120	{
3121	b = tree->a_edges[br];
3122
3123	old_pars = tree->c_pars;
3124	Spr_Subtree(b,b->left,tree);
3125	new_pars = tree->c_pars;
3126
3127	pars_diff = new_pars - old_pars;
3128	if(pars_diff > max_pars_diff) max_pars_diff = pars_diff;
3129
3130	old_pars = tree->c_pars;
3131	Spr_Subtree(b,b->rght,tree);
3132	new_pars = tree->c_pars;
3133
3134	pars_diff = new_pars - old_pars;
3135	if(pars_diff > max_pars_diff) max_pars_diff = pars_diff;
3136	}
3137
3138	/* tree->mod->s_opt->pars_thresh = MAX(5,max_pars_diff); */
3139
3140	return 1;
3141	}
3142
3143	/*********************************************************/
3144
3145	void Spr_Subtree(t_edge b, t_node link, t_tree *tree)
3146	{
3147	int i;
3148	int n_moves_pars, n_moves, curr_pars, min_pars, best_move;
3149	t_spr *best_pars_move;
3150	t_edge target, residual;
3151
3152	best_move = -1;
3153	tree->n_moves = 0;
3154	curr_pars = tree->c_pars;
3155
3156	MIXT_Set_Pars_Thresh(tree);
3157
3158	if((link != b->left) && (link != b->rght))
3159	{
3160	PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
3161	Exit("\n");
3162	}
3163	else
3164	{
3165	/* printf("\n. -1"); fflush(NULL); */
3166	/* if(!Check_Lk_At_Given_Edge(NO,tree)) Exit("\n"); */
3167
3168	if(!link->tax) Test_All_Spr_Targets(b,link,tree);
3169
3170	/* printf("\n. 0"); fflush(NULL); */
3171	/* if(!Check_Lk_At_Given_Edge(NO,tree)) Exit("\n"); */
3172
3173	if(tree->n_moves)
3174	{
3175	n_moves_pars = 0;
3176	n_moves = 0;
3177
3178	For(i,tree->n_moves)
3179	if(curr_pars - tree->spr_list[i]->pars >= -tree->mod->s_opt->pars_thresh)
3180	n_moves_pars++;
3181
3182
3183	For(i,tree->n_moves)
3184	{
3185	n_moves++;
3186	/* if(n_moves > 15) break; */
3187	if(n_moves > 5) break;
3188	if(tree->spr_list[i]->lnL < tree->best_lnL - 2. * tree->mod->s_opt->max_delta_lnL_spr) break;
3189	}
3190
3191	if(!tree->mod->s_opt->spr_lnL) n_moves = n_moves_pars;
3192
3193	if(!tree->io->fp_in_constraint_tree) n_moves = MAX(1,n_moves);
3194	n_moves = MIN(n_moves,2*tree->n_otu-3);
3195
3196	if(tree->mod->s_opt->spr_pars)
3197	{
3198	if(tree->io->fp_in_constraint_tree)
3199	{
3200	PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
3201	Exit("\n");
3202	}
3203
3204	min_pars = 1E+8;
3205	best_pars_move = NULL;
3206
3207
3208	For(i,n_moves)
3209	if(tree->spr_list[i]->pars < min_pars)
3210	{
3211	best_pars_move = tree->spr_list[i];
3212	min_pars = tree->spr_list[i]->pars;
3213	}
3214
3215	if(best_pars_move->pars < tree->best_pars)
3216	{
3217	Prune_Subtree(best_pars_move->n_link,best_pars_move->n_opp_to_link,&target,&residual,tree);
3218	Graft_Subtree(best_pars_move->b_target,best_pars_move->n_link,residual,tree);
3219	Set_Both_Sides(YES,tree);
3220	Pars(NULL,tree);
3221	tree->best_pars = tree->c_pars;
3222	if(tree->best_pars != best_pars_move->pars)
3223	{
3224	PhyML_Printf("\n== best_pars = %d move_pars = %d",tree->best_pars,best_pars_move->pars);
3225	PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
3226	Exit("\n");
3227	}
3228	tree->n_improvements++;
3229	}
3230	else
3231	Pars(NULL,tree);
3232	}
3233	else
3234	{
3235	/* tree->both_sides = YES; */
3236	/* Lk(NULL,tree); */
3237	/* tree->both_sides = NO; */
3238
3239	/* printf("\n. 1"); fflush(NULL); */
3240	/* if(!Check_Lk_At_Given_Edge(YES,tree)) Exit("\n"); */
3241
3242	best_move = Evaluate_List_Of_Regraft_Pos_Triple(tree->spr_list,n_moves,tree);
3243
3244	/* printf("\n. 2"); fflush(NULL); */
3245	/* if(!Check_Lk_At_Given_Edge(YES,tree)) Exit("\n"); */
3246
3247	if((best_move > -1) && (tree->spr_list[best_move]->lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move))
3248	/* if((best_move > -1) && (tree->spr_list[best_move]->lnL > tree->best_lnL - tree->mod->s_opt->max_delta_lnL_spr)) */
3249	{
3250	Try_One_Spr_Move_Triple(tree->spr_list[best_move],tree);
3251	}
3252	else
3253	{
3254	Pars(NULL,tree);
3255	/* tree->both_sides = YES; */
3256	/* Lk(tree); */
3257	/* tree->both_sides = NO; */
3258	}
3259
3260	/* printf("\n. 3"); fflush(NULL); */
3261	/* if(!Check_Lk_At_Given_Edge(tree)) Exit("\n"); */
3262
3263	}
3264	}
3265	Reset_Spr_List(tree);
3266	}
3267	}
3268
3269
3270	/*********************************************************/
3271
3272	int Test_All_Spr_Targets(t_edge b_pulled, t_node n_link, t_tree *tree)
3273	{
3274	t_node n_opp_to_link,n_v1,*n_v2;
3275	t_edge b_target,b_residual;
3276	int i,dir1,dir2;
3277	phydbl init_len_v1, init_len_v2, *init_len_pulled;
3278	int best_found;
3279	phydbl init_lnL;
3280
3281	if(tree->mixt_tree != NULL)
3282	{
3283	PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
3284	Exit("\n");
3285	}
3286
3287
3288
3289	init_lnL = tree->c_lnL;
3290	b_target = b_residual = NULL;
3291	n_opp_to_link = (n_link == b_pulled->rght)?(b_pulled->left):(b_pulled->rght);
3292
3293	init_len_pulled = MIXT_Get_Lengths_Of_This_Edge(b_pulled,tree);
3294
3295	dir1 = dir2 = -1;
3296	For(i,3)
3297	if(n_link->v[i] != n_opp_to_link)
3298	{
3299	if(dir1<0) dir1 = i;
3300	else dir2 = i;
3301	}
3302
3303	if(n_link->v[dir1]->num < n_link->v[dir2]->num)
3304	{
3305	n_v1 = n_link->v[dir1];
3306	n_v2 = n_link->v[dir2];
3307	init_len_v1 = MIXT_Get_Lengths_Of_This_Edge(n_link->b[dir1],tree);
3308	init_len_v2 = MIXT_Get_Lengths_Of_This_Edge(n_link->b[dir2],tree);
3309	}
3310	else
3311	{
3312	n_v1 = n_link->v[dir2];
3313	n_v2 = n_link->v[dir1];
3314	init_len_v1 = MIXT_Get_Lengths_Of_This_Edge(n_link->b[dir2],tree);
3315	init_len_v2 = MIXT_Get_Lengths_Of_This_Edge(n_link->b[dir1],tree);
3316	}
3317
3318	if(!(n_v1->tax && n_v2->tax)) /! Pruning is meaningless otherwise /
3319	{
3320
3321	Prune_Subtree(n_link,n_opp_to_link,&b_target,&b_residual,tree);
3322
3323	if(tree->mod->s_opt->spr_lnL)
3324	{
3325	Fast_Br_Len(b_target,tree,NO);
3326	/* Update_PMat_At_Given_Edge(b_target,tree); */
3327	}
3328
3329
3330	best_found = 0;
3331	tree->depth_curr_path = 0;
3332	tree->curr_path[0] = b_target->left;
3333	Test_One_Spr_Target_Recur(b_target->rght,
3334	b_target->left,
3335	b_pulled,n_link,b_residual,&best_found,tree);
3336
3337
3338	tree->depth_curr_path = 0;
3339	tree->curr_path[0] = b_target->rght;
3340	Test_One_Spr_Target_Recur(b_target->left,
3341	b_target->rght,
3342	b_pulled,n_link,b_residual,&best_found,tree);
3343
3344
3345	Graft_Subtree(b_target,n_link,b_residual,tree);
3346
3347
3348	if((n_link->v[dir1] != n_v1) \|\| (n_link->v[dir2] != n_v2))
3349	PhyML_Printf("\n. Warning: -- SWITCH NEEDED -- ! \n");
3350
3351	/* n_link->b[dir1]->l->v = init_len_v1; */
3352	/* n_link->b[dir2]->l->v = init_len_v2; */
3353	/* b_pulled->l->v = init_len_pulled; */
3354
3355	MIXT_Set_Lengths_Of_This_Edge(init_len_v1,n_link->b[dir1],tree);
3356	MIXT_Set_Lengths_Of_This_Edge(init_len_v2,n_link->b[dir2],tree);
3357	MIXT_Set_Lengths_Of_This_Edge(init_len_pulled,b_pulled,tree);
3358
3359	Update_PMat_At_Given_Edge(n_link->b[dir1],tree);
3360	Update_PMat_At_Given_Edge(n_link->b[dir2],tree);
3361	Update_PMat_At_Given_Edge(b_pulled,tree);
3362
3363	/* if(tree->mod->s_opt->spr_lnL) */
3364	/* { */
3365	/* I don't understand why this is required when spr_lnL = NO, but it is... */
3366	MIXT_Set_Alias_Subpatt(YES,tree);
3367	Update_P_Lk(tree,b_pulled, n_link);
3368	Update_P_Lk(tree,b_target, n_link);
3369	Update_P_Lk(tree,b_residual,n_link);
3370	MIXT_Set_Alias_Subpatt(NO,tree);
3371	/* } */
3372	/* else */
3373	/* { */
3374	Update_P_Pars(tree,b_pulled, n_link);
3375	Update_P_Pars(tree,b_target, n_link);
3376	Update_P_Pars(tree,b_residual,n_link);
3377	/* } */
3378
3379	For(i,3)
3380	if(n_link->v[i] != n_opp_to_link)
3381	{
3382	if(tree->mod->s_opt->spr_lnL)
3383	{
3384	MIXT_Set_Alias_Subpatt(YES,tree);
3385	Pre_Order_Lk(n_link,n_link->v[i],tree);
3386	MIXT_Set_Alias_Subpatt(NO,tree);
3387	}
3388	else
3389	Pre_Order_Pars(n_link,n_link->v[i],tree);
3390	}
3391
3392	/* printf("\n. 000000"); fflush(NULL); */
3393	/* if(!Check_Lk_At_Given_Edge(NO,tree)) Exit("\n"); */
3394	}
3395
3396	tree->c_lnL = init_lnL;
3397
3398	Free(init_len_v1);
3399	Free(init_len_v2);
3400	Free(init_len_pulled);
3401
3402	return 0;
3403
3404	}
3405
3406	/*********************************************************/
3407
3408	void Test_One_Spr_Target_Recur(t_node a, t_node d, t_edge pulled, t_node link, t_edge residual, int best_found, t_tree *tree)
3409	{
3410	int i;
3411
3412	if(*best_found) return;
3413
3414	if(d->tax) return;
3415	else
3416	{
3417	phydbl move_lnL;
3418
3419	For(i,3)
3420	{
3421	if(d->v[i] != a)
3422	{
3423
3424	if(tree->mod->s_opt->spr_lnL)
3425	{
3426	MIXT_Set_Alias_Subpatt(YES,tree);
3427	Update_P_Lk(tree,d->b[i],d);
3428	MIXT_Set_Alias_Subpatt(NO,tree);
3429	}
3430	else
3431	Update_P_Pars(tree,d->b[i],d);
3432
3433	tree->depth_curr_path++;
3434	tree->curr_path[tree->depth_curr_path] = d->v[i];
3435
3436	if((tree->depth_curr_path <= tree->mod->s_opt->max_depth_path) &&
3437	(tree->depth_curr_path >= tree->mod->s_opt->min_depth_path))
3438	{
3439	move_lnL = Test_One_Spr_Target(d->b[i],pulled,link,residual,tree);
3440	if(move_lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move)
3441	{
3442	*best_found = 1;
3443	}
3444	}
3445
3446	if(tree->depth_curr_path < tree->mod->s_opt->max_depth_path)
3447	Test_One_Spr_Target_Recur(d,d->v[i],pulled,link,residual,best_found,tree);
3448
3449	tree->depth_curr_path--;
3450	}
3451	}
3452	}
3453	}
3454
3455	/*********************************************************/
3456
3457	phydbl Test_One_Spr_Target(t_edge b_target, t_edge b_arrow, t_node n_link, t_edge b_residual, t_tree *tree)
3458	{
3459	phydbl init_target_len, init_arrow_len, *init_residual_len;
3460	int i,dir_v0,dir_v1,dir_v2;
3461	phydbl l0,l1,*l2;
3462	t_node v1, v2;
3463	phydbl init_lnL, move_lnL;
3464	int init_pars;
3465	t_tree *orig_tree;
3466	t_edge orig_target, orig_arrow;
3467	t_node *orig_link;
3468	t_spr orig_move,move;
3469
3470	if(tree->mixt_tree != NULL)
3471	{
3472	PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
3473	Exit("\n");
3474	}
3475
3476	tree->n_moves++;
3477
3478	move_lnL = UNLIKELY;
3479	init_lnL = tree->c_lnL;
3480	init_pars = tree->c_pars;
3481
3482	Graft_Subtree(b_target,n_link,b_residual,tree);
3483
3484	/* init_target_len = b_target->l->v; */
3485	/* init_arrow_len = b_arrow->l->v; */
3486	/* init_residual_len = b_residual->l->v; */
3487
3488	init_target_len = MIXT_Get_Lengths_Of_This_Edge(b_target,tree);
3489	init_arrow_len = MIXT_Get_Lengths_Of_This_Edge(b_arrow,tree);
3490	init_residual_len = MIXT_Get_Lengths_Of_This_Edge(b_residual,tree);
3491
3492	if(tree->mod->s_opt->spr_lnL)
3493	{
3494	MIXT_Set_Alias_Subpatt(YES,tree);
3495	/* move_lnL = Triple_Dist(n_link,tree,1); */
3496	Update_PMat_At_Given_Edge(b_target,tree);
3497	Update_PMat_At_Given_Edge(b_arrow,tree);
3498	Update_P_Lk(tree,b_residual,n_link);
3499	move_lnL = Lk(b_residual,tree);
3500	MIXT_Set_Alias_Subpatt(NO,tree);
3501	}
3502	else
3503	{
3504	Update_P_Pars(tree,b_residual,n_link);
3505	Pars(b_residual,tree);
3506	}
3507
3508	v1 = (b_residual->left == n_link)?(b_residual->rght):(b_residual->left);
3509	v2 = (b_target->left == n_link)?(b_target->rght):(b_target->left);
3510	dir_v1 = dir_v2 = dir_v0 = -1;
3511	For(i,3)
3512	{
3513	if(n_link->v[i] == v1) dir_v1 = i;
3514	else if(n_link->v[i] == v2) dir_v2 = i;
3515	else dir_v0 = i;
3516	}
3517
3518	l0 = MIXT_Get_Lengths_Of_This_Edge(n_link->b[dir_v0],tree);
3519	if(n_link->v[dir_v1]->num > n_link->v[dir_v2]->num)
3520	{
3521	l1 = MIXT_Get_Lengths_Of_This_Edge(n_link->b[dir_v2],tree);
3522	l2 = MIXT_Get_Lengths_Of_This_Edge(n_link->b[dir_v1],tree);
3523	}
3524	else
3525	{
3526	l1 = MIXT_Get_Lengths_Of_This_Edge(n_link->b[dir_v1],tree);
3527	l2 = MIXT_Get_Lengths_Of_This_Edge(n_link->b[dir_v2],tree);
3528	}
3529
3530	move = tree->spr_list[tree->size_spr_list];
3531	For(i,tree->depth_curr_path+1) move->path[i] = tree->curr_path[i];
3532	move->depth_path = tree->depth_curr_path;
3533	move->pars = tree->c_pars;
3534	move->lnL = tree->c_lnL;
3535
3536	orig_target = b_target;
3537	orig_link = n_link;
3538	orig_arrow = b_arrow;
3539	orig_tree = tree;
3540	orig_move = move;
3541	i = 0;
3542	do
3543	{
3544	move->l0 = l0[i];
3545	move->l1 = l1[i];
3546	move->l2 = l2[i];
3547
3548	if(tree->mod->gamma_mgf_bl == YES)
3549	{
3550	move->v0 = l0[i+1];
3551	move->v1 = l1[i+1];
3552	move->v2 = l2[i+1];
3553	}
3554
3555	move->b_target = b_target;
3556	move->n_link = n_link;
3557	move->b_opp_to_link = b_arrow;
3558	move->dist = b_target->topo_dist_btw_edges;
3559	move->n_opp_to_link = (n_link==b_arrow->left)?(b_arrow->rght):(b_arrow->left);
3560
3561	if(tree->next)
3562	{
3563	tree = tree->next;
3564	b_target = b_target->next;
3565	b_arrow = b_arrow->next;
3566	n_link = n_link->next;
3567	move = move->next;
3568	}
3569	else
3570	{
3571	tree = tree->next;
3572	b_target = b_target->next;
3573	b_arrow = b_arrow->next;
3574	n_link = n_link->next;
3575	move = move->next;
3576	}
3577
3578	i+=2;
3579	}
3580	while(tree);
3581
3582	b_target = orig_target;
3583	n_link = orig_link;
3584	b_arrow = orig_arrow;
3585	tree = orig_tree;
3586	move = orig_move;
3587
3588	Include_One_Spr_To_List_Of_Spr(move,tree);
3589
3590	/* b_target->l->v = init_target_len; */
3591	/* b_arrow->l->v = init_arrow_len; */
3592	/* b_residual->l->v = init_residual_len; */
3593
3594	MIXT_Set_Lengths_Of_This_Edge(init_target_len,b_target,tree);
3595	MIXT_Set_Lengths_Of_This_Edge(init_arrow_len,b_arrow,tree);
3596	MIXT_Set_Lengths_Of_This_Edge(init_residual_len,b_residual,tree);
3597
3598	Prune_Subtree(n_link,
3599	(n_link==b_arrow->left)?(b_arrow->rght):(b_arrow->left),
3600	&b_target,
3601	&b_residual,
3602	tree);
3603
3604	if(tree->mod->s_opt->spr_lnL) Update_PMat_At_Given_Edge(b_target,tree);
3605
3606	tree->c_lnL = init_lnL;
3607	tree->c_pars = init_pars;
3608
3609	Free(l0);
3610	Free(l1);
3611	Free(l2);
3612
3613	Free(init_target_len);
3614	Free(init_arrow_len);
3615	Free(init_residual_len);
3616
3617	return move_lnL;
3618	}
3619
3620	/*********************************************************/
3621
3622	void Speed_Spr_Loop(t_tree *tree)
3623	{
3624	phydbl lk_old;
3625
3626	tree->best_pars = 1E+8;
3627	tree->mod->s_opt->spr_lnL = 0;
3628	tree->mod->s_opt->spr_pars = 0;
3629	tree->mod->s_opt->quickdirty = 0;
3630
3631	if((tree->mod->s_opt->print) && (!tree->io->quiet)) PhyML_Printf("\n\n. Maximizing likelihood (using SPR moves)...\n");
3632
3633	SPR_Shuffle(tree);
3634
3635	Optimiz_All_Free_Param(tree,(tree->io->quiet)?(0):(tree->mod->s_opt->print));
3636	tree->best_lnL = tree->c_lnL;
3637
3638	/*****************************/
3639	lk_old = UNLIKELY;
3640	tree->mod->s_opt->max_depth_path = 2*tree->n_otu-3;
3641	tree->mod->s_opt->max_delta_lnL_spr = (tree->io->datatype == NT)?(10.):(0.);
3642	/* tree->mod->s_opt->max_delta_lnL_spr = (tree->io->datatype == NT)?(50.):(0.); */
3643	/* tree->mod->s_opt->max_depth_path = 5; */
3644	tree->mod->s_opt->spr_lnL = NO;
3645	do
3646	{
3647	lk_old = tree->c_lnL;
3648	Speed_Spr(tree,1);
3649	if(tree->n_improvements) Optimiz_All_Free_Param(tree,(tree->io->quiet)?(0):(tree->mod->s_opt->print));
3650	if((tree->io->datatype == NT) && ((tree->max_spr_depth < 4 && tree->n_improvements < 10) \|\| (FABS(lk_old-tree->c_lnL) < 1.))) break;
3651	if((tree->io->datatype == AA) && ((tree->max_spr_depth < 4 && tree->n_improvements < 1) \|\| (FABS(lk_old-tree->c_lnL) < 1.))) break;
3652	/* if((tree->io->datatype == NT) & (tree->n_improvements < 10) \|\| (FABS(lk_old-tree->c_lnL) < 1.)) break; */
3653	/* if((tree->io->datatype == AA) & (tree->n_improvements < 1) \|\| (FABS(lk_old-tree->c_lnL) < 1.)) break; */
3654	}
3655	while(1);
3656	/*****************************/
3657
3658
3659	/*****************************/
3660	if(tree->io->datatype == NT)
3661	{
3662	lk_old = UNLIKELY;
3663	tree->mod->s_opt->max_delta_lnL_spr = 20.;
3664	tree->mod->s_opt->max_depth_path = 10;
3665	tree->mod->s_opt->spr_lnL = YES;
3666	do
3667	{
3668	lk_old = tree->c_lnL;
3669	Speed_Spr(tree,1);
3670	if(tree->n_improvements) Optimiz_All_Free_Param(tree,(tree->io->quiet)?(0):(tree->mod->s_opt->print));
3671	if((!tree->n_improvements) \|\| (FABS(lk_old-tree->c_lnL) < 1.)) break;
3672	}
3673	while(1);
3674	}
3675	/*****************************/
3676
3677
3678	/*****************************/
3679	lk_old = UNLIKELY;
3680	do
3681	{
3682	lk_old = tree->c_lnL;
3683	if(!Simu(tree,10)) break;
3684	Optimiz_All_Free_Param(tree,(tree->io->quiet)?(0):(tree->mod->s_opt->print));
3685	}
3686	while(FABS(lk_old - tree->c_lnL) > tree->mod->s_opt->min_diff_lk_local);
3687	/*****************************/
3688
3689
3690	/*****************************/
3691	do
3692	{
3693	Round_Optimize(tree,tree->data,ROUND_MAX);
3694	if(!Check_NNI_Five_Branches(tree)) break;
3695	}while(1);
3696	/*****************************/
3697
3698	/* if((tree->mod->s_opt->print) && (!tree->io->quiet)) PhyML_Printf("\n"); */
3699
3700	}
3701
3702	//////////////////////////////////////////////////////////////
3703	//////////////////////////////////////////////////////////////
3704
3705	void Speed_Spr(t_tree *tree, int max_cycles)
3706	{
3707	int step,old_pars;
3708	phydbl old_lnL;
3709
3710	if(tree->lock_topo == YES)
3711	{
3712	PhyML_Printf("\n== The tree topology is locked.");
3713	PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
3714	Exit("\n");
3715	}
3716
3717	Set_Both_Sides(YES,tree);
3718	Pars(NULL,tree);
3719	Lk(NULL,tree);
3720	Record_Br_Len(tree);
3721
3722	tree->mod->s_opt->deepest_path = 0;
3723	tree->best_pars = tree->c_pars;
3724	tree->best_lnL = tree->c_lnL;
3725	old_lnL = tree->c_lnL;
3726	old_pars = tree->c_pars;
3727	step = 0;
3728	do
3729	{
3730	++step;
3731
3732	old_lnL = tree->c_lnL;
3733	old_pars = tree->c_pars;
3734
3735	tree->max_spr_depth = 0;
3736	tree->n_improvements = 0;
3737	tree->perform_spr_right_away = 1;
3738	Spr(UNLIKELY,tree);
3739
3740	if(!tree->mod->s_opt->spr_pars)
3741	{
3742	/* Optimise branch lengths */
3743	Optimize_Br_Len_Serie(tree);
3744	/* Update partial likelihoods */
3745	Set_Both_Sides(YES,tree);
3746	Lk(NULL,tree);
3747
3748	/* Print log-likelihood and parsimony scores */
3749	if((tree->mod->s_opt->print) && (!tree->io->quiet)) Print_Lk(tree,"[Branch lengths ]");
3750	}
3751	else
3752	{
3753	if((tree->mod->s_opt->print) && (!tree->io->quiet)) Print_Pars(tree);
3754	}
3755
3756	Pars(NULL,tree);
3757
3758	if(tree->io->print_trace)
3759	{
3760	char *s = Write_Tree(tree,NO);
3761	PhyML_Fprintf(tree->io->fp_out_trace,"[%f]%s\n",tree->c_lnL,s); fflush(tree->io->fp_out_trace);
3762	if((tree->io->print_site_lnl) && (!tree->mod->s_opt->spr_pars)) { Print_Site_Lk(tree,tree->io->fp_out_lk); } fflush(tree->io->fp_out_lk);
3763	Free(s);
3764	}
3765
3766	/* Record the current best log-likelihood and parsimony */
3767	tree->best_lnL = tree->c_lnL;
3768	tree->best_pars = tree->c_pars;
3769
3770	if(!tree->mod->s_opt->spr_pars)
3771	{
3772	if(tree->c_lnL < old_lnL-tree->mod->s_opt->min_diff_lk_local)
3773	{
3774	PhyML_Printf("\n== old_lnL = %f c_lnL = %f",old_lnL,tree->c_lnL);
3775	PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
3776	Exit("");
3777	}
3778	}
3779	else
3780	{
3781	if(tree->c_pars > old_pars)
3782	{
3783	PhyML_Printf("\n== old_pars = %d c_pars = %d",old_pars,tree->c_pars);
3784	PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
3785	Exit("");
3786	}
3787	}
3788
3789	/* Record the current best branch lengths */
3790	Record_Br_Len(tree);
3791
3792	/* Exit if no improvements after complete optimization */
3793	if(step+1 > max_cycles) break;
3794	if((!tree->mod->s_opt->spr_pars) && (FABS(old_lnL-tree->c_lnL) < tree->mod->s_opt->min_diff_lk_global)) break;
3795	/* if(( tree->mod->s_opt->spr_pars) && (FABS(old_pars-tree->c_pars) < 1.)) break; */
3796	if(!tree->n_improvements) break;
3797
3798	}while(1);
3799	}
3800
3801	/*********************************************************/
3802
3803	int Evaluate_List_Of_Regraft_Pos_Triple(t_spr *spr_list, int list_size, t_tree tree)
3804	{
3805	t_spr move,orig_move;
3806	t_edge init_target, b_residual;
3807	int i,j,best_move,n;
3808	int dir_v0, dir_v1, dir_v2;
3809	phydbl *recorded_l;
3810	phydbl best_lnL,init_lnL;
3811	int recorded;
3812	t_tree *orig_tree;
3813
3814	if(tree->mixt_tree != NULL)
3815	{
3816	PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
3817	Exit("\n");
3818	}
3819
3820	best_lnL = UNLIKELY;
3821	init_target = b_residual = NULL;
3822	best_move = -1;
3823	init_lnL = tree->c_lnL;
3824	recorded_l = NULL;
3825
3826	if(!list_size && !tree->io->fp_in_constraint_tree)
3827	{
3828	PhyML_Printf("\n== List size is 0 !");
3829	PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
3830	Exit("\n");
3831	}
3832
3833	recorded = NO;
3834	For(i,list_size)
3835	{
3836	move = spr_list[i];
3837
3838	if(!move)
3839	{
3840	PhyML_Printf("\n== move is NULL\n");
3841	PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
3842	Exit("\n");
3843	}
3844
3845	if(move->b_target)
3846	{
3847	/* Record t_edge lengths */
3848	Record_Br_Len(tree);
3849
3850	/* Prune subtree */
3851	Prune_Subtree(move->n_link,move->n_opp_to_link,&init_target,&b_residual,tree);
3852
3853	if(recorded == NO)
3854	{
3855	/*! Rough optimisation of the branch length at prune site
3856	* We only need to perform this optimisation for the first
3857	* element of spr_list because the pruned subtree is the
3858	* same across all the elements of spr_list. It would not
3859	* be true in the general case
3860	*/
3861	recorded = YES;
3862
3863	Fast_Br_Len(init_target,tree,NO);
3864
3865	/! Record branch length at prune site /
3866	recorded_l = MIXT_Get_Lengths_Of_This_Edge(init_target,tree);
3867
3868	/! Make sure recorded_l has been updated beforehand /
3869	orig_move = move;
3870	orig_tree = tree;
3871	n = 0;
3872	do
3873	{
3874	move->init_target_l = recorded_l[n];
3875	move->init_target_v = recorded_l[n+1];
3876	n+=2;
3877
3878
3879	if(tree->next)
3880	{
3881	move = move->next;
3882	tree = tree->next;
3883	}
3884	else
3885	{
3886	move = move->next;
3887	tree = tree->next;
3888	}
3889	}
3890	while(tree);
3891	move = orig_move;
3892	tree = orig_tree;
3893	}
3894	else
3895	{
3896	MIXT_Set_Lengths_Of_This_Edge(recorded_l,init_target,tree);
3897
3898	orig_move = move;
3899	orig_tree = tree;
3900	n = 0;
3901	do
3902	{
3903	move->init_target_l = recorded_l[n];
3904	move->init_target_v = recorded_l[n+1];
3905	n+=2;
3906
3907	if(tree->next)
3908	{
3909	move = move->next;
3910	tree = tree->next;
3911	}
3912	else
3913	{
3914	move = move->next;
3915	tree = tree->next;
3916	}
3917	}
3918	while(tree);
3919	move = orig_move;
3920	tree = orig_tree;
3921	}
3922
3923	/* Update the change proba matrix at prune position */
3924	Update_PMat_At_Given_Edge(init_target,tree);
3925
3926	/* Update conditional likelihoods along the path from the prune to
3927	the regraft position */
3928	MIXT_Set_Alias_Subpatt(YES,tree);
3929	Update_P_Lk_Along_A_Path(move->path,move->depth_path+1,tree);
3930	MIXT_Set_Alias_Subpatt(NO,tree);
3931
3932	/* Regraft subtree */
3933	Graft_Subtree(move->b_target,move->n_link,b_residual,tree);
3934
3935	dir_v1 = dir_v2 = dir_v0 = -1;
3936	For(j,3)
3937	{
3938	if(move->n_link->v[j] == move->n_opp_to_link) dir_v0 = j;
3939	else if(dir_v1 < 0) dir_v1 = j;
3940	else dir_v2 = j;
3941	}
3942
3943	orig_tree = tree;
3944	orig_move = move;
3945	do
3946	{
3947	move->n_link->b[dir_v0]->l->v = move->l0;
3948	move->n_link->b[dir_v0]->l_var->v = move->v0;
3949
3950	if(move->n_link->v[dir_v1]->num > move->n_link->v[dir_v2]->num)
3951	{
3952	move->n_link->b[dir_v2]->l->v = move->l1;
3953	move->n_link->b[dir_v1]->l->v = move->l2;
3954
3955	if(tree->io->mod->gamma_mgf_bl == YES)
3956	{
3957	move->n_link->b[dir_v2]->l_var->v = move->v1;
3958	move->n_link->b[dir_v1]->l_var->v = move->v2;
3959	}
3960	}
3961	else
3962	{
3963	move->n_link->b[dir_v1]->l->v = move->l1;
3964	move->n_link->b[dir_v2]->l->v = move->l2;
3965
3966	if(tree->io->mod->gamma_mgf_bl == YES)
3967	{
3968	move->n_link->b[dir_v1]->l_var->v = move->v1;
3969	move->n_link->b[dir_v2]->l_var->v = move->v2;
3970	}
3971	}
3972
3973	if(tree->next)
3974	{
3975	move = move->next;
3976	tree = tree->next;
3977	}
3978	else
3979	{
3980	move = move->next;
3981	tree = tree->next;
3982	}
3983
3984	}
3985	while(tree);
3986	move = orig_move;
3987	tree = orig_tree;
3988
3989	MIXT_Set_Alias_Subpatt(YES,tree);
3990	move->lnL = Triple_Dist(move->n_link,tree,YES);
3991	MIXT_Set_Alias_Subpatt(NO,tree);
3992
3993	if((move->lnL < best_lnL) && (move->lnL > best_lnL - tree->mod->s_opt->max_delta_lnL_spr))
3994	{
3995	/* Estimate the three t_edge lengths at the regraft site */
3996	move->lnL = Triple_Dist(move->n_link,tree,NO);
3997	}
3998
3999	/* Record updated branch lengths for this move */
4000	orig_move = move;
4001	orig_tree = tree;
4002	do
4003	{
4004	move->l0 = move->n_link->b[dir_v0]->l->v;
4005	move->v0 = move->n_link->b[dir_v0]->l_var->v;
4006
4007	if(move->n_link->v[dir_v1]->num > move->n_link->v[dir_v2]->num)
4008	{
4009	move->l1 = move->n_link->b[dir_v2]->l->v;
4010	move->l2 = move->n_link->b[dir_v1]->l->v;
4011
4012	if(tree->io->mod->gamma_mgf_bl == YES)
4013	{
4014	move->v1 = move->n_link->b[dir_v2]->l_var->v;
4015	move->v2 = move->n_link->b[dir_v1]->l_var->v;
4016	}
4017	}
4018	else
4019	{
4020	move->l1 = move->n_link->b[dir_v1]->l->v;
4021	move->l2 = move->n_link->b[dir_v2]->l->v;
4022
4023	if(tree->io->mod->gamma_mgf_bl == YES)
4024	{
4025	move->v1 = move->n_link->b[dir_v1]->l_var->v;
4026	move->v2 = move->n_link->b[dir_v2]->l_var->v;
4027	}
4028	}
4029
4030	/* printf("\n. %f %f %f %f",move->l0,move->l1,move->l2,move->lnL); */
4031	/* For(j,2tree->n_otu-3) /
4032	/* printf("\n== %d %f %f", */
4033	/* j, */
4034	/* tree->a_edges[j]->gamma_prior_mean, */
4035	/* tree->a_edges[j]->l_var->v); */
4036
4037	if(tree->next)
4038	{
4039	move = move->next;
4040	tree = tree->next;
4041	}
4042	else
4043	{
4044	move = move->next;
4045	tree = tree->next;
4046	}
4047	}
4048	while(tree);
4049	move = orig_move;
4050	tree = orig_tree;
4051
4052	if(move->lnL > best_lnL + tree->mod->s_opt->min_diff_lk_move)
4053	{
4054	best_lnL = move->lnL;
4055	best_move = i;
4056	}
4057
4058	/* Regraft the subtree at its original position */
4059	Prune_Subtree(move->n_link,
4060	move->n_opp_to_link,
4061	&move->b_target,
4062	&b_residual,
4063	tree);
4064
4065	Graft_Subtree(init_target,
4066	move->n_link,
4067	b_residual,
4068	tree);
4069
4070	/* Restore branch lengths */
4071	Restore_Br_Len(tree);
4072
4073	/* Update relevant change proba matrices */
4074	/* Update_PMat_At_Given_Edge(move->n_link->b[0],tree); */
4075	/* Update_PMat_At_Given_Edge(move->n_link->b[1],tree); */
4076	/* Update_PMat_At_Given_Edge(move->n_link->b[2],tree); */
4077	Update_PMat_At_Given_Edge(move->b_target,tree);
4078
4079	/* Update_P_Lk(tree,move->n_link->b[0],move->n_link); */
4080	/* Update_P_Lk(tree,move->n_link->b[1],move->n_link); */
4081	/* Update_P_Lk(tree,move->n_link->b[2],move->n_link); */
4082
4083	/* Update conditional likelihoods along the path from the prune to
4084	the regraft position */
4085	/* Update_P_Lk_Along_A_Path(move->path,move->depth_path+1,tree); */
4086
4087
4088	tree->c_lnL = init_lnL;
4089	}
4090
4091	/* Bail out as soon as you've found a true improvement */
4092	if(move->lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move) break;
4093	}
4094
4095	/* PhyML_Printf("\n. [ %4d/%4d ]",i,list_size); */
4096	/* PhyML_Printf("\n. max_improv = %f",max_improv); */
4097
4098
4099	MIXT_Set_Alias_Subpatt(YES,tree);
4100	For(i,list_size)
4101	{
4102	move = spr_list[i];
4103	if(move->b_target)
4104	{
4105
4106	For(j,3) Update_PMat_At_Given_Edge(move->n_link->b[j],tree);
4107	For(j,3) Update_P_Lk(tree,move->n_link->b[j],move->n_link);
4108
4109	/* TO DO : we don't need to update all these partial likelihoods here.
4110	Would need to record only those that were along the paths examined
4111	above */
4112
4113	For(j,3)
4114	if(move->n_link->v[j] != move->n_opp_to_link)
4115	Pre_Order_Lk(move->n_link,move->n_link->v[j],tree);
4116
4117	break;
4118	}
4119	}
4120
4121	MIXT_Set_Alias_Subpatt(NO,tree);
4122
4123	#ifdef DEBUG
4124	if(best_move < 0 && list_size > 0)
4125	{
4126	PhyML_Printf("\n\n. Best_move < 0 !");
4127
4128	PhyML_Printf("\n. List size = %d",list_size);
4129	For(i,list_size)
4130	{
4131	move = spr_list[i];
4132	PhyML_Printf("\n. %p %p",move,move->b_target);
4133	}
4134
4135	PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
4136	Exit("\n");
4137	}
4138	#endif
4139
4140	Free(recorded_l);
4141
4142	return best_move;
4143	}
4144
4145	/*********************************************************/
4146
4147	int Try_One_Spr_Move_Triple(t_spr move, t_tree tree)
4148	{
4149	t_edge init_target, b_residual;
4150	int j;
4151	int dir_v0, dir_v1, dir_v2;
4152	int accept;
4153	t_edge *orig_b;
4154	t_spr *orig_move;
4155	t_tree *orig_tree;
4156
4157	if(tree->mixt_tree != NULL)
4158	{
4159	PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
4160	Exit("\n");
4161	}
4162
4163	Record_Br_Len(tree);
4164
4165	Prune_Subtree(move->n_link,
4166	move->n_opp_to_link,
4167	&init_target,
4168	&b_residual,
4169	tree);
4170
4171	orig_tree = tree;
4172	orig_b = init_target;
4173	orig_move = move;
4174	do
4175	{
4176	init_target->l->v = move->init_target_l;
4177	init_target->l_var->v = move->init_target_v;
4178
4179	if(tree->next)
4180	{
4181	init_target = init_target->next;
4182	move = move->next;
4183	tree = tree->next;
4184	}
4185	else
4186	{
4187	init_target = init_target->next;
4188	move = move->next;
4189	tree = tree->next;
4190	}
4191	}
4192	while(tree);
4193	init_target = orig_b;
4194	move = orig_move;
4195	tree = orig_tree;
4196
4197	Graft_Subtree(move->b_target,move->n_link,b_residual,tree);
4198
4199	dir_v1 = dir_v2 = dir_v0 = -1;
4200	For(j,3)
4201	{
4202	if(move->n_link->v[j] == move->n_opp_to_link) dir_v0 = j;
4203	else if(dir_v1 < 0) dir_v1 = j;
4204	else dir_v2 = j;
4205	}
4206
4207	orig_move = move;
4208	orig_tree = tree;
4209	do
4210	{
4211	move->n_link->b[dir_v0]->l->v = move->l0;
4212	move->n_link->b[dir_v0]->l_var->v = move->v0;
4213
4214	if(move->n_link->v[dir_v1]->num > move->n_link->v[dir_v2]->num)
4215	{
4216	move->n_link->b[dir_v2]->l->v = move->l1;
4217	move->n_link->b[dir_v1]->l->v = move->l2;
4218
4219	if(tree->io->mod->gamma_mgf_bl == YES)
4220	{
4221	move->n_link->b[dir_v2]->l_var->v = move->v1;
4222	move->n_link->b[dir_v1]->l_var->v = move->v2;
4223	}
4224	}
4225	else
4226	{
4227	move->n_link->b[dir_v1]->l->v = move->l1;
4228	move->n_link->b[dir_v2]->l->v = move->l2;
4229
4230	if(tree->io->mod->gamma_mgf_bl == YES)
4231	{
4232	move->n_link->b[dir_v1]->l_var->v = move->v1;
4233	move->n_link->b[dir_v2]->l_var->v = move->v2;
4234	}
4235	}
4236
4237	if(tree->next)
4238	{
4239	move = move->next;
4240	tree = tree->next;
4241	}
4242	else
4243	{
4244	move = move->next;
4245	tree = tree->next;
4246	}
4247
4248	}
4249	while(tree);
4250	move = orig_move;
4251	tree = orig_tree;
4252
4253	accept = YES;
4254	if(!Check_Topo_Constraints(tree,tree->io->cstr_tree)) accept = NO;
4255
4256	if(accept == YES) /* Apply the move */
4257	{
4258	time(&(tree->t_current));
4259	Pars(NULL,tree);
4260	Set_Both_Sides(YES,tree);
4261	MIXT_Set_Alias_Subpatt(YES,tree);
4262	Lk(NULL,tree);
4263	MIXT_Set_Alias_Subpatt(NO,tree);
4264
4265	if(FABS(tree->c_lnL - move->lnL) > tree->mod->s_opt->min_diff_lk_move)
4266	{
4267	int i;
4268	PhyML_Printf("\n== c_lnL = %f move_lnL = %f", tree->c_lnL,move->lnL);
4269	printf("\n== l0=%f l1=%f l2=%f v0=%f v1=%f v2=%f",move->l0,move->l1,move->l2,move->v0,move->v1,move->v2);
4270	For(i,2*tree->n_otu-3)
4271	printf("\n== %d %f %f",
4272	i,
4273	tree->a_edges[i]->l->v,
4274	tree->a_edges[i]->l_var->v);
4275
4276	/* printf("\n. %f %f %f %f", */
4277	/* tree->next->c_lnL, */
4278	/* tree->next->next->c_lnL, */
4279	/* tree->next->next->c_lnL, */
4280	/* tree->next->next->next->c_lnL); */
4281	PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
4282	Exit("\n");
4283	}
4284
4285	//
4286	/* if(!(tree->n_improvements % tree->mod->s_opt->br_len_in_spr)) */
4287	/* { */
4288	/* tree->mod->s_opt->brent_it_max = 10; */
4289	/* Optimize_Br_Len_Serie(tree->a_nodes[0], */
4290	/* tree->a_nodes[0]->v[0], */
4291	/* tree->a_nodes[0]->b[0], */
4292	/* tree, */
4293	/* tree->data); */
4294	/* tree->mod->s_opt->brent_it_max = 500; */
4295
4296	/* tree->both_sides = 1; */
4297	/* Lk(tree); */
4298	/* } */
4299	//
4300
4301
4302	if((tree->mod->s_opt->print) && (!tree->io->quiet))
4303	{
4304	Print_Lk_And_Pars(tree);
4305	PhyML_Printf(" [depth=%5d]",move->depth_path); fflush(NULL);
4306	}
4307
4308	if(move->depth_path > tree->max_spr_depth) tree->max_spr_depth = move->depth_path;
4309
4310	tree->n_improvements++;
4311	tree->best_lnL = tree->c_lnL;
4312	Record_Br_Len(tree);
4313
4314	if(move->depth_path > tree->mod->s_opt->deepest_path)
4315	tree->mod->s_opt->deepest_path = move->depth_path;
4316
4317	return 1;
4318	}
4319
4320	Prune_Subtree(move->n_link,
4321	move->n_opp_to_link,
4322	&move->b_target,
4323	&b_residual,
4324	tree);
4325
4326	Graft_Subtree(init_target,
4327	move->n_link,
4328	b_residual,
4329	tree);
4330
4331	Restore_Br_Len(tree);
4332
4333	Set_Both_Sides(YES,tree);
4334	MIXT_Set_Alias_Subpatt(YES,tree);
4335	Lk(NULL,tree);
4336	MIXT_Set_Alias_Subpatt(NO,tree);
4337	Pars(NULL,tree);
4338	return 0;
4339	}
4340
4341	/*********************************************************/
4342
4343	int Try_One_Spr_Move_Full(t_spr move, t_tree tree)
4344	{
4345	t_edge init_target, b_residual;
4346
4347	Record_Br_Len(tree);
4348
4349	Prune_Subtree(move->n_link,
4350	move->n_opp_to_link,
4351	&init_target,
4352	&b_residual,
4353	tree);
4354
4355	Graft_Subtree(move->b_target,move->n_link,b_residual,tree);
4356
4357	MIXT_Set_Alias_Subpatt(YES,tree);
4358	Set_Both_Sides(YES,tree);
4359	Lk(NULL,tree);
4360	MIXT_Set_Alias_Subpatt(NO,tree);
4361
4362	Optimize_Br_Len_Serie(tree);
4363
4364	Set_Both_Sides(YES,tree);
4365	Lk(NULL,tree);
4366
4367	if(tree->c_lnL > tree->best_lnL + tree->mod->s_opt->min_diff_lk_move)
4368	{
4369	Pars(NULL,tree);
4370	if((tree->mod->s_opt->print) && (!tree->io->quiet)) Print_Lk(tree,"[Topology ]");
4371	tree->n_improvements++;
4372	tree->best_lnL = tree->c_lnL;
4373	Record_Br_Len(tree);
4374	return 1;
4375	}
4376	else
4377	{
4378	Prune_Subtree(move->n_link,
4379	move->n_opp_to_link,
4380	&move->b_target,
4381	&b_residual,
4382	tree);
4383
4384	Graft_Subtree(init_target,
4385	move->n_link,
4386	b_residual,
4387	tree);
4388
4389	Restore_Br_Len(tree);
4390	Set_Both_Sides(YES,tree);
4391
4392	MIXT_Set_Alias_Subpatt(YES,tree);
4393	Lk(NULL,tree);
4394	MIXT_Set_Alias_Subpatt(NO,tree);
4395	Pars(NULL,tree);
4396	return 0;
4397	}
4398
4399	return -1;
4400	}
4401
4402	/*********************************************************/
4403
4404	void Include_One_Spr_To_List_Of_Spr(t_spr move, t_tree tree)
4405	{
4406	int i;
4407	t_spr buff_spr,orig_move, orig_move_list, move_list;
4408	t_tree *orig_tree;
4409
4410	if(tree->mixt_tree != NULL)
4411	{
4412	PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
4413	Exit("\n");
4414	}
4415
4416	if((( tree->mod->s_opt->spr_lnL) && (move->lnL > tree->spr_list[tree->size_spr_list-1]->lnL)) \|\|
4417	((!tree->mod->s_opt->spr_lnL) && (move->pars <= tree->spr_list[tree->size_spr_list-1]->pars)))
4418	{
4419	move_list = tree->spr_list[tree->size_spr_list-1];
4420
4421	move_list->depth_path = move->depth_path;
4422	move_list->pars = move->pars;
4423	move_list->lnL = move->lnL;
4424	move_list->dist = move->dist;
4425
4426	orig_move = move;
4427	orig_move_list = move_list;
4428	orig_tree = tree;
4429	do
4430	{
4431	move_list->l0 = move->l0;
4432	move_list->l1 = move->l1;
4433	move_list->l2 = move->l2;
4434	move_list->v0 = move->v0;
4435	move_list->v1 = move->v1;
4436	move_list->v2 = move->v2;
4437	move_list->b_target = move->b_target;
4438	move_list->n_link = move->n_link;
4439	move_list->n_opp_to_link = move->n_opp_to_link;
4440	move_list->b_opp_to_link = move->b_opp_to_link;
4441
4442
4443	if(tree->next)
4444	{
4445	move = move->next;
4446	move_list = move_list->next;
4447	tree = tree->next;
4448	}
4449	else
4450	{
4451	move = move->next;
4452	move_list = move_list->next;
4453	tree = tree->next;
4454	}
4455	}
4456	while(tree);
4457	move = orig_move;
4458	move_list = orig_move_list;
4459	tree = orig_tree;
4460
4461	For(i,move_list->depth_path+1) move_list->path[i] = move->path[i];
4462
4463	for(i=tree->size_spr_list-1;i>0;i--)
4464	{
4465	if((( tree->mod->s_opt->spr_lnL) && (tree->spr_list[i]->lnL > tree->spr_list[i-1]->lnL)) \|\|
4466	((!tree->mod->s_opt->spr_lnL) && (tree->spr_list[i]->pars <= tree->spr_list[i-1]->pars)))
4467	{
4468
4469	orig_tree = tree;
4470	do
4471	{
4472	buff_spr = tree->spr_list[i-1];
4473	tree->spr_list[i-1] = tree->spr_list[i];
4474	tree->spr_list[i] = buff_spr;
4475
4476	if(tree->next) tree = tree->next;
4477	else tree = tree->next;
4478	}
4479	while(tree);
4480	tree = orig_tree;
4481
4482	}
4483	else break;
4484	}
4485	}
4486	}
4487
4488	/*********************************************************/
4489
4490	void Random_Spr(int n_moves, t_tree *tree)
4491	{
4492	int i;
4493	int br_pulled, br_target;
4494	t_spr *spr_struct;
4495	t_edge target, residual;
4496
4497	spr_struct = Make_One_Spr(tree);
4498	Init_One_Spr(spr_struct);
4499	target = residual = NULL;
4500	For(i,n_moves)
4501	{
4502	br_pulled = (int)((phydbl)rand()/RAND_MAX * (2*tree->n_otu-3-1));
4503	do
4504	{
4505	br_target = (int)((phydbl)rand()/RAND_MAX * (2*tree->n_otu-3-1));
4506	}while(br_target == br_pulled);
4507
4508	spr_struct->n_link = tree->a_edges[br_pulled]->left;
4509	spr_struct->n_opp_to_link = tree->a_edges[br_pulled]->rght;
4510	spr_struct->b_opp_to_link = tree->a_edges[br_pulled];
4511	spr_struct->b_target = tree->a_edges[br_target];
4512	spr_struct->b_init_target = NULL;
4513
4514	if(!Check_Spr_Move_Validity(spr_struct,tree))
4515	{
4516	spr_struct->n_link = tree->a_edges[br_pulled]->rght;
4517	spr_struct->n_opp_to_link = tree->a_edges[br_pulled]->left;
4518	}
4519
4520	#ifdef DEBUG
4521	if(!Check_Spr_Move_Validity(spr_struct,tree))
4522	{
4523	Warn_And_Exit("\n. Could not find a valid move...\n");
4524	}
4525	#endif
4526
4527	Prune_Subtree(spr_struct->n_link,
4528	spr_struct->n_opp_to_link,
4529	&target,
4530	&residual,
4531	tree);
4532
4533	Graft_Subtree(spr_struct->b_target,
4534	spr_struct->n_link,
4535	residual,tree);
4536	}
4537	Free(spr_struct);
4538	}
4539
4540	/*********************************************************/
4541
4542	void Reset_Spr_List(t_tree *tree)
4543	{
4544	int i;
4545
4546	For(i,tree->size_spr_list)
4547	{
4548	tree->spr_list[i]->depth_path = 0;
4549	tree->spr_list[i]->pars = MAX_PARS;
4550	tree->spr_list[i]->lnL = UNLIKELY;
4551	tree->spr_list[i]->n_link = NULL;
4552	tree->spr_list[i]->n_opp_to_link = NULL;
4553	tree->spr_list[i]->b_target = NULL;
4554	}
4555	}
4556
4557	/*********************************************************/
4558
4559	int Check_Spr_Move_Validity(t_spr this_spr_move, t_tree tree)
4560	{
4561	int match;
4562
4563	match = 0;
4564	Found_In_Subtree(this_spr_move->n_link,
4565	this_spr_move->n_opp_to_link,
4566	this_spr_move->b_target->left,
4567	&match,
4568	tree);
4569
4570	if(match) return 0;
4571	else return 1;
4572	}
4573
4574	/*********************************************************/
4575
4576	void Spr_Pars(t_tree *tree)
4577	{
4578
4579	PhyML_Printf("\n. Minimizing parsimony...\n");
4580
4581	tree->best_pars = 1E+8;
4582	tree->best_lnL = UNLIKELY;
4583	tree->mod->s_opt->spr_lnL = 0;
4584
4585	tree->mod->s_opt->spr_pars = 1;
4586	do
4587	{
4588	Speed_Spr(tree,1);
4589	}while(tree->n_improvements);
4590	tree->mod->s_opt->spr_pars = 0;
4591
4592	PhyML_Printf("\n");
4593	}
4594
4595	//////////////////////////////////////////////////////////////
4596	//////////////////////////////////////////////////////////////
4597
4598
4599	void SPR_Shuffle(t_tree *mixt_tree)
4600	{
4601	phydbl lk_old;
4602	int *orig_catg,n;
4603	t_tree tree,*tree_list;
4604
4605	if(mixt_tree->mod->s_opt->print) PhyML_Printf("\n\n. Refining the tree...\n");
4606
4607	/! Get the number of classes in each mixture /
4608	orig_catg = MIXT_Get_Number_Of_Classes_In_All_Mixtures(mixt_tree);
4609
4610
4611	/*! Set the number of rate classes to (at most) 2.
4612	! Propagate this to every mixture tree in the analysis
4613	*/
4614	tree = mixt_tree;
4615	n = 0;
4616	do
4617	{
4618	tree->mod->ras->n_catg = MIN(2,orig_catg[n]);
4619	if(tree->mod->ras->invar == YES) tree->mod->ras->n_catg--;
4620	tree = tree->next_mixt;
4621	n++;
4622	}
4623	while(tree);
4624
4625
4626	/*! Make sure the number of trees in each mixture is at most 2
4627	*/
4628	tree_list = MIXT_Record_All_Mixtures(mixt_tree);
4629	MIXT_Break_All_Mixtures(orig_catg,mixt_tree);
4630
4631	Set_Both_Sides(YES,mixt_tree);
4632	Lk(NULL,mixt_tree);
4633
4634
4635	/* mixt_tree->mod->s_opt->print = YES; */
4636	mixt_tree->best_pars = 1E+8;
4637	mixt_tree->mod->s_opt->spr_pars = NO;
4638	mixt_tree->mod->s_opt->quickdirty = NO;
4639	mixt_tree->best_lnL = mixt_tree->c_lnL;
4640	mixt_tree->mod->s_opt->max_delta_lnL_spr = 0.;
4641	mixt_tree->mod->s_opt->max_depth_path = 2*mixt_tree->n_otu-3;
4642	mixt_tree->mod->s_opt->spr_lnL = NO;
4643
4644	do
4645	{
4646	Set_Both_Sides(YES,mixt_tree);
4647	Lk(NULL,mixt_tree);
4648	Pars(NULL,mixt_tree);
4649	Record_Br_Len(mixt_tree);
4650
4651	mixt_tree->n_improvements = 0;
4652	mixt_tree->max_spr_depth = 0;
4653	mixt_tree->best_pars = mixt_tree->c_pars;
4654	mixt_tree->best_lnL = mixt_tree->c_lnL;
4655
4656	lk_old = mixt_tree->c_lnL;
4657	Spr(UNLIKELY,mixt_tree);
4658
4659	Optimiz_All_Free_Param(mixt_tree,(mixt_tree->io->quiet)?(0):(mixt_tree->mod->s_opt->print));
4660
4661	Optimize_Br_Len_Serie(mixt_tree);
4662
4663	if(mixt_tree->n_improvements < 20 \|\| mixt_tree->max_spr_depth < 5 \|\|
4664	(FABS(lk_old-mixt_tree->c_lnL) < 1.)) break;
4665	}
4666	while(1);
4667
4668
4669	if(mixt_tree->mod->s_opt->print && (!mixt_tree->io->quiet))
4670	PhyML_Printf("\n\n. End of refining stage...\n. The log-likelihood might now decrease and then increase again...\n");
4671
4672
4673	/*! Go back to the original data structure, with potentially more
4674	! than 2 trees per mixture
4675	*/
4676	MIXT_Reconnect_All_Mixtures(tree_list,mixt_tree);
4677	Free(tree_list);
4678
4679	/*! Set the number of rate classes to their original values
4680	*/
4681	tree = mixt_tree;
4682	n = 0;
4683	do
4684	{
4685	tree->mod->ras->n_catg = orig_catg[n];
4686	if(tree->mod->ras->invar == YES) tree->mod->ras->n_catg--;
4687	tree = tree->next_mixt;
4688	n++;
4689	}
4690	while(tree);
4691
4692	Free(orig_catg);
4693
4694	/*! Only the first two trees for each mixture have been modified so
4695	! far -> need to update the other trees by copying the modified trees
4696	! onto them.
4697	*/
4698	tree = mixt_tree;
4699	do
4700	{
4701	if(tree != mixt_tree) Copy_Tree(mixt_tree,tree);
4702	tree = tree->next;
4703	}
4704	while(tree);
4705	}
4706
4707
4708
4709
4710
4711
4712
4713
4714
4715
4716	/*
4717	** EOF: spr.c
4718	*/

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

Context Navigation

source: trunk/GDE/PHYML20130708/phyml/src/spr.c

Download in other formats: