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 #include "make.h"
14
15 //////////////////////////////////////////////////////////////
16
Make_Tree_For_Lk(t_tree * tree)17 void Make_Tree_For_Lk(t_tree *tree)
18 {
19 int i;
20 calign *cdata;
21
22 cdata = tree->data;
23 assert(cdata);
24
25
26 tree->c_lnL_sorted = (phydbl *)mCalloc(tree->n_pattern,sizeof(phydbl));
27 tree->cur_site_lk = (phydbl *)mCalloc(tree->n_pattern,sizeof(phydbl));
28 tree->old_site_lk = (phydbl *)mCalloc(tree->n_pattern,sizeof(phydbl));
29 tree->site_lk_cat = (phydbl *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(phydbl));
30 tree->unscaled_site_lk_cat = (phydbl *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->n_pattern,sizeof(phydbl));
31 tree->fact_sum_scale = (int *)mCalloc(tree->n_pattern,sizeof(int));
32
33
34 #if (defined(__AVX__) || defined(__SSE3__))
35 #ifndef WIN32
36 if(posix_memalign((void **)&tree->eigen_lr_left,BYTE_ALIGN,(size_t)MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
37 if(posix_memalign((void **)&tree->eigen_lr_rght,BYTE_ALIGN,(size_t)MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
38 if(posix_memalign((void **)&tree->dot_prod,BYTE_ALIGN,(size_t)tree->n_pattern*tree->mod->ns*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
39 if(posix_memalign((void **)&tree->expl,BYTE_ALIGN,(size_t)3*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
40 #else
41 tree->eigen_lr_left = _aligned_malloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
42 tree->eigen_lr_rght = _aligned_malloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
43 tree->dot_prod = _aligned_malloc(tree->n_pattern*tree->mod->ns*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*sizeof(phydbl),BYTE_ALIGN);
44 tree->expl = _aligned_malloc(3*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
45 #endif
46 #else
47 tree->eigen_lr_left = (phydbl *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns,sizeof(phydbl));
48 tree->eigen_lr_rght = (phydbl *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns,sizeof(phydbl));
49 tree->dot_prod = (phydbl *)mCalloc(tree->n_pattern*tree->mod->ns*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(phydbl));
50 tree->expl = (phydbl *)mCalloc(3*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns,sizeof(phydbl));
51 #endif
52
53
54 tree->log_lks_aLRT = (phydbl **)mCalloc(3,sizeof(phydbl *));
55 for(i=0;i<3;i++) tree->log_lks_aLRT[i] = (phydbl *)mCalloc(tree->data->init_len,sizeof(phydbl));
56
57 for(i=0;i<2*tree->n_otu-1;++i) Make_Edge_NNI(tree->a_edges[i]);
58
59 Make_Extra_Edge_Lk(tree);
60
61 if(tree->is_mixt_tree == NO)
62 {
63 for(i=0;i<2*tree->n_otu-1;++i) Make_Edge_Lk(tree->a_edges[i],tree);
64 for(i=0;i<2*tree->n_otu-2;++i) Make_Node_Lk(tree->a_nodes[i]);
65 for(i=0;i<2*tree->n_otu-1;++i) Make_Edge_Loc(tree->a_edges[i],tree);
66
67 Init_Partial_Lk_Tips_Double(tree);
68 Init_Partial_Lk_Loc(tree);
69
70 if(tree->n_root != NULL)
71 {
72 Free_Edge_Lk_Rght(tree->n_root->b[1]);
73 Free_Edge_Lk_Rght(tree->n_root->b[2]);
74 Free_Edge_Loc_Rght(tree->n_root->b[1]);
75 Free_Edge_Loc_Rght(tree->n_root->b[2]);
76
77 tree->n_root->b[1]->p_lk_rght = tree->e_root->p_lk_left;
78 tree->n_root->b[2]->p_lk_rght = tree->e_root->p_lk_rght;
79
80 tree->n_root->b[1]->p_lk_tip_r = tree->e_root->p_lk_tip_l;
81 tree->n_root->b[2]->p_lk_tip_r = tree->e_root->p_lk_tip_r;
82
83 tree->n_root->b[1]->div_post_pred_rght = tree->e_root->div_post_pred_rght;
84 tree->n_root->b[2]->div_post_pred_rght = tree->e_root->div_post_pred_left;
85
86 tree->n_root->b[1]->sum_scale_rght = tree->e_root->sum_scale_rght;
87 tree->n_root->b[2]->sum_scale_rght = tree->e_root->sum_scale_left;
88
89 tree->n_root->b[1]->sum_scale_rght_cat = tree->e_root->sum_scale_rght_cat;
90 tree->n_root->b[2]->sum_scale_rght_cat = tree->e_root->sum_scale_left_cat;
91
92 tree->n_root->b[1]->patt_id_rght = tree->e_root->patt_id_rght;
93 tree->n_root->b[2]->patt_id_rght = tree->e_root->patt_id_left;
94 }
95 }
96 }
97
98 //////////////////////////////////////////////////////////////
99 //////////////////////////////////////////////////////////////
100
Make_Tree_For_Pars(t_tree * tree)101 void Make_Tree_For_Pars(t_tree *tree)
102 {
103 int i;
104 calign *cdata;
105
106 cdata = tree->data;
107 assert(cdata);
108
109 assert(tree->mod);
110
111 tree->site_pars = (int *)mCalloc(tree->n_pattern,sizeof(int));
112 tree->step_mat = (int *)mCalloc(tree->mod->ns * tree->mod->ns,sizeof(int));
113
114 for(i=0;i<2*tree->n_otu-1;++i) Make_Edge_Pars(tree->a_edges[i],tree);
115
116 Init_Ui_Tips(tree);
117 Init_Partial_Pars_Tips(tree); /* Must be called after Init_Ui_Tips is called */
118
119 if(tree->n_root)
120 {
121 Free_Edge_Pars_Rght(tree->a_edges[2*tree->n_otu-3]);
122 Free_Edge_Pars_Rght(tree->a_edges[2*tree->n_otu-2]);
123 }
124
125 Get_Step_Mat(tree);
126 }
127
128 //////////////////////////////////////////////////////////////
129 //////////////////////////////////////////////////////////////
130
Make_All_Edges_Lk(t_node * a,t_node * d,t_tree * tree)131 void Make_All_Edges_Lk(t_node *a, t_node *d, t_tree *tree)
132 {
133 int i;
134
135 for(i=0;i<3;i++)
136 if((a->v[i]) && (a->v[i] == d))
137 Make_Edge_Lk(a->b[i],tree);
138
139 if(d->tax) return;
140 else
141 {
142 for(i=0;i<3;i++)
143 {
144 if(d->v[i] != a && d->b[i] != tree->e_root)
145 Make_All_Edges_Lk(d,d->v[i],tree);
146 }
147 }
148 }
149
150 //////////////////////////////////////////////////////////////
151 //////////////////////////////////////////////////////////////
152
Make_Edge_Light(t_node * a,t_node * d,int num)153 t_edge *Make_Edge_Light(t_node *a, t_node *d, int num)
154 {
155 t_edge *b;
156
157 b = (t_edge *)mCalloc(1,sizeof(t_edge));
158
159 b->l = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
160 Init_Scalar_Dbl(b->l);
161
162 b->l_old = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
163 Init_Scalar_Dbl(b->l_old);
164
165 b->l_var = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
166 Init_Scalar_Dbl(b->l_var);
167
168 b->l_var_old = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
169 Init_Scalar_Dbl(b->l_var_old);
170
171 Init_Edge_Light(b,num);
172
173 if(a && b)
174 {
175 b->left = a;
176 b->rght = d;
177 if(a->tax) {b->rght = a; b->left = d;} /* root */
178 /* a tip is necessary on the right side of the t_edge */
179
180 (b->left == a)?
181 (Set_Edge_Dirs(b,a,d,NULL)):
182 (Set_Edge_Dirs(b,d,a,NULL));
183
184 assert(b->l_r > -1);
185 assert(b->r_l > -1);
186
187 b->l_old->v = b->l->v;
188 }
189 else
190 {
191 b->left = NULL;
192 b->rght = NULL;
193 }
194
195 return b;
196
197 }
198
199 //////////////////////////////////////////////////////////////
200 //////////////////////////////////////////////////////////////
201
Make_Edge_Pars(t_edge * b,t_tree * tree)202 void Make_Edge_Pars(t_edge *b, t_tree *tree)
203 {
204 assert(b);
205 Make_Edge_Pars_Left(b,tree);
206 Make_Edge_Pars_Rght(b,tree);
207 }
208
209 //////////////////////////////////////////////////////////////
210 //////////////////////////////////////////////////////////////
211
Make_Edge_Pars_Left(t_edge * b,t_tree * tree)212 void Make_Edge_Pars_Left(t_edge *b, t_tree *tree)
213 {
214 b->pars_l = (int *)mCalloc(tree->data->crunch_len,sizeof(int));
215 b->ui_l = (int *)mCalloc(tree->data->crunch_len,sizeof(int));
216 b->p_pars_l = (int *)mCalloc(tree->data->crunch_len*tree->mod->ns,sizeof(int ));
217 b->n_diff_states_l = (int *)mCalloc(tree->mod->ns,sizeof(int ));
218 }
219
220 //////////////////////////////////////////////////////////////
221 //////////////////////////////////////////////////////////////
222
Make_Edge_Pars_Rght(t_edge * b,t_tree * tree)223 void Make_Edge_Pars_Rght(t_edge *b, t_tree *tree)
224 {
225 b->pars_r = (int *)mCalloc(tree->data->crunch_len,sizeof(int));
226 b->ui_r = (int *)mCalloc(tree->data->crunch_len,sizeof(int));
227 b->p_pars_r = (int *)mCalloc(tree->data->crunch_len*tree->mod->ns,sizeof(int ));
228 b->n_diff_states_r = (int *)mCalloc(tree->mod->ns,sizeof(int ));
229 }
230
231 //////////////////////////////////////////////////////////////
232 //////////////////////////////////////////////////////////////
233
Make_Edge_Lk(t_edge * b,t_tree * tree)234 void Make_Edge_Lk(t_edge *b, t_tree *tree)
235 {
236 if(tree->is_mixt_tree)
237 {
238 PhyML_Printf("\n== Err. in file %s at line %d (function '%s') \n",__FILE__,__LINE__,__FUNCTION__);
239 Warn_And_Exit("");
240 }
241
242 b->l_old->v = b->l->v;
243
244 #if (defined(__AVX__) || defined(__SSE3__))
245 #ifndef WIN32
246 if(posix_memalign((void *)&b->Pij_rr,BYTE_ALIGN,(size_t)tree->mod->ras->n_catg*tree->mod->ns*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
247 if(posix_memalign((void *)&b->tPij_rr,BYTE_ALIGN,(size_t)tree->mod->ras->n_catg*tree->mod->ns*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
248 #else
249 b->Pij_rr = _aligned_malloc(tree->mod->ras->n_catg*tree->mod->ns*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
250 b->tPij_rr = _aligned_malloc(tree->mod->ras->n_catg*tree->mod->ns*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
251 #endif
252 #else
253 b->Pij_rr = (phydbl *)mCalloc(tree->mod->ras->n_catg*tree->mod->ns*tree->mod->ns,sizeof(phydbl));
254 b->tPij_rr = (phydbl *)mCalloc(tree->mod->ras->n_catg*tree->mod->ns*tree->mod->ns,sizeof(phydbl));
255 #endif
256
257 Make_Edge_Lk_Left(b,tree);
258 Make_Edge_Lk_Rght(b,tree);
259 }
260
261 //////////////////////////////////////////////////////////////
262 //////////////////////////////////////////////////////////////
263
Make_Edge_Lk_Left(t_edge * b,t_tree * tree)264 void Make_Edge_Lk_Left(t_edge *b, t_tree *tree)
265 {
266 int ns = tree->mod->ns;
267
268 b->div_post_pred_left = (short int *)mCalloc(ns,sizeof(short int));
269
270 b->sum_scale_left_cat = (int *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int));
271
272 if(b->left && !b->left->tax)
273 b->sum_scale_left = (int *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int));
274 else
275 b->sum_scale_left = NULL;
276
277 if(b->left)
278 {
279 if((!b->left->tax) || (tree->mod->s_opt->greedy))
280 {
281 #if (defined(__AVX__) || defined(__SSE3__))
282 #ifndef WIN32
283 if(posix_memalign((void **)&b->p_lk_left,BYTE_ALIGN,(size_t)tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
284 #else
285 b->p_lk_left = _aligned_malloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
286 #endif
287 #else
288 b->p_lk_left = (phydbl *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns,sizeof(phydbl));
289 #endif
290 b->p_lk_tip_l = NULL;
291 }
292 else if(b->left->tax)
293 {
294 b->p_lk_left = NULL;
295
296 #if (defined(__AVX__) || defined(__SSE3__))
297 #ifndef WIN32
298 if(posix_memalign((void **)&b->p_lk_tip_l,BYTE_ALIGN,(size_t)tree->data->crunch_len*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
299 #else
300 b->p_lk_tip_l = _aligned_malloc(tree->data->crunch_len*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
301 #endif
302 #else
303 b->p_lk_tip_l = (phydbl *)mCalloc(tree->data->crunch_len*tree->mod->ns,sizeof(phydbl));
304 #endif
305 }
306 }
307 else
308 {
309 b->p_lk_left = NULL;
310 b->p_lk_tip_l = NULL;
311 }
312
313 if(b->num >= 2*tree->n_otu-3)
314 {
315 b->sum_scale_left = (int *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int));
316
317 #if (defined(__AVX__) || defined(__SSE3__))
318 #ifndef WIN32
319 if(posix_memalign((void **)&b->p_lk_left,BYTE_ALIGN,(size_t)tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
320 #else
321 b->p_lk_left = _aligned_malloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
322 #endif
323 #else
324 b->p_lk_left = (phydbl *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns,sizeof(phydbl));
325 #endif
326 }
327
328 b->patt_id_left = (int *)mCalloc(tree->data->crunch_len,sizeof(int));
329 }
330
331 //////////////////////////////////////////////////////////////
332 //////////////////////////////////////////////////////////////
333
Make_Edge_Lk_Rght(t_edge * b,t_tree * tree)334 void Make_Edge_Lk_Rght(t_edge *b, t_tree *tree)
335 {
336 int ns = tree->mod->ns;
337
338 b->div_post_pred_rght = (short int *)mCalloc(ns,sizeof(short int));
339
340 b->sum_scale_rght_cat = (int *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int));
341
342 if(b->rght && !b->rght->tax)
343 b->sum_scale_rght = (int *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int));
344 else
345 b->sum_scale_rght = NULL;
346
347
348 if(b->rght)
349 {
350 if((!b->rght->tax) || (tree->mod->s_opt->greedy))
351 {
352 #if (defined(__AVX__) || defined(__SSE3__))
353 #ifndef WIN32
354 if(posix_memalign((void **)&b->p_lk_rght,BYTE_ALIGN,(size_t)tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
355 #else
356 b->p_lk_rght = _aligned_malloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
357 #endif
358 #else
359 b->p_lk_rght = (phydbl *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns,sizeof(phydbl));
360 #endif
361 b->p_lk_tip_r = NULL;
362 }
363 else if(b->rght->tax)
364 {
365 #if (defined(__AVX__) || defined(__SSE3__))
366 #ifndef WIN32
367 if(posix_memalign((void **)&b->p_lk_tip_r,BYTE_ALIGN,(size_t)tree->data->crunch_len*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
368 #else
369 b->p_lk_tip_r = _aligned_malloc(tree->data->crunch_len*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
370 #endif
371 #else
372 b->p_lk_tip_r = (phydbl *)mCalloc(tree->data->crunch_len*tree->mod->ns,sizeof(phydbl));
373 #endif
374 b->p_lk_rght = NULL;
375 }
376 }
377 else
378 {
379 b->p_lk_rght = NULL;
380 b->p_lk_tip_r = NULL;
381 }
382
383 if(b->num >= 2*tree->n_otu-3)
384 {
385 b->sum_scale_rght = (int *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int));
386 #if (defined(__AVX__) || defined(__SSE3__))
387 #ifndef WIN32
388 if(posix_memalign((void **)&b->p_lk_rght,BYTE_ALIGN,(size_t)tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
389 #else
390 b->p_lk_rght = _aligned_malloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
391 #endif
392
393
394 #else
395 b->p_lk_rght = (phydbl *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns,sizeof(phydbl));
396 #endif
397 }
398
399 b->patt_id_rght = (int *)mCalloc(tree->data->crunch_len,sizeof(int));
400 }
401
402 //////////////////////////////////////////////////////////////
403 //////////////////////////////////////////////////////////////
404
Make_Extra_Edge_Lk(t_tree * tree)405 void Make_Extra_Edge_Lk(t_tree *tree)
406 {
407 int ns = tree->mod->ns;
408
409 tree->div_post_pred_extra_0 = (short int *)mCalloc(ns,sizeof(short int));
410
411 tree->sum_scale_cat_extra_0 = (int *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int));
412
413 tree->sum_scale_extra_0 = (int *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int));
414
415 #if (defined(__AVX__) || defined(__SSE3__))
416 #ifndef WIN32
417 if(posix_memalign((void **)&tree->p_lk_extra_0,BYTE_ALIGN,(size_t)tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
418 #else
419 tree->p_lk_extra_0 = _aligned_malloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
420 #endif
421 #else
422 tree->p_lk_extra_0 = (phydbl *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns,sizeof(phydbl));
423 #endif
424
425
426 #if (defined(__AVX__) || defined(__SSE3__))
427 #ifndef WIN32
428 if(posix_memalign((void **)&tree->p_lk_tip_extra_0,BYTE_ALIGN,(size_t)tree->data->crunch_len*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
429 #else
430 tree->p_lk_tip_extra_0 = _aligned_malloc(tree->data->crunch_len*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
431 #endif
432 #else
433 tree->p_lk_tip_extra_0 = (phydbl *)mCalloc(tree->data->crunch_len*tree->mod->ns,sizeof(phydbl));
434 #endif
435
436 tree->patt_id_extra_0 = (int *)mCalloc(tree->data->crunch_len,sizeof(int));
437
438 tree->div_post_pred_extra_1 = (short int *)mCalloc(ns,sizeof(short int));
439
440 tree->sum_scale_cat_extra_1 = (int *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int));
441
442 tree->sum_scale_extra_1 = (int *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int));
443
444 #if (defined(__AVX__) || defined(__SSE3__))
445 #ifndef WIN32
446 if(posix_memalign((void **)&tree->p_lk_extra_1,BYTE_ALIGN,(size_t)tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
447 #else
448 tree->p_lk_extra_1 = _aligned_malloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
449 #endif
450 #else
451 tree->p_lk_extra_1 = (phydbl *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)*tree->mod->ns,sizeof(phydbl));
452 #endif
453
454
455 #if (defined(__AVX__) || defined(__SSE3__))
456 #ifndef WIN32
457 if(posix_memalign((void **)&tree->p_lk_tip_extra_1,BYTE_ALIGN,(size_t)tree->data->crunch_len*tree->mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
458 #else
459 tree->p_lk_tip_extra_1 = _aligned_malloc(tree->data->crunch_len*tree->mod->ns*sizeof(phydbl),BYTE_ALIGN);
460 #endif
461 #else
462 tree->p_lk_tip_extra_1 = (phydbl *)mCalloc(tree->data->crunch_len*tree->mod->ns,sizeof(phydbl));
463 #endif
464
465 tree->patt_id_extra_1 = (int *)mCalloc(tree->data->crunch_len,sizeof(int));
466
467 }
468
469 //////////////////////////////////////////////////////////////
470 //////////////////////////////////////////////////////////////
471
Make_Edge_Loc(t_edge * b,t_tree * tree)472 void Make_Edge_Loc(t_edge *b, t_tree *tree)
473 {
474 Make_Edge_Loc_Left(b,tree);
475 Make_Edge_Loc_Rght(b,tree);
476 }
477
478 //////////////////////////////////////////////////////////////
479 //////////////////////////////////////////////////////////////
480
Make_Edge_Loc_Rght(t_edge * b,t_tree * tree)481 void Make_Edge_Loc_Rght(t_edge *b, t_tree *tree)
482 {
483 b->p_lk_loc_rght = (int *)mCalloc(tree->data->crunch_len,sizeof(int));
484 }
485
486 //////////////////////////////////////////////////////////////
487 //////////////////////////////////////////////////////////////
488
Make_Edge_Loc_Left(t_edge * b,t_tree * tree)489 void Make_Edge_Loc_Left(t_edge *b, t_tree *tree)
490 {
491 b->p_lk_loc_left = (int *)mCalloc(tree->data->crunch_len,sizeof(int));
492 }
493
494 //////////////////////////////////////////////////////////////
495 //////////////////////////////////////////////////////////////
496
Make_Edge_NNI(t_edge * b)497 void Make_Edge_NNI(t_edge *b)
498 {
499 b->nni = Make_NNI();
500 b->nni->b = b;
501 b->nni->left = b->left;
502 b->nni->rght = b->rght;
503 }
504
505 //////////////////////////////////////////////////////////////
506 //////////////////////////////////////////////////////////////
507
Make_NNI()508 t_nni *Make_NNI()
509 {
510 t_nni *a_nni;
511
512 a_nni = (t_nni *)mCalloc(1,sizeof(t_nni));
513
514 Init_NNI(a_nni);
515 return a_nni;
516 }
517
518 //////////////////////////////////////////////////////////////
519 //////////////////////////////////////////////////////////////
520
Make_Node_Light(int num)521 t_node *Make_Node_Light(int num)
522 {
523 t_node *n;
524
525 n = (t_node *)mCalloc(1,sizeof(t_node));
526 n->v = (t_node **)mCalloc(3,sizeof(t_node *));
527 n->b = (t_edge **)mCalloc(3,sizeof(t_edge *));
528 n->score = (phydbl *)mCalloc(3,sizeof(phydbl));
529 n->s_ingrp = (int *)mCalloc(3,sizeof(int));
530 n->s_outgrp = (int *)mCalloc(3,sizeof(int));
531 n->cal = (t_cal **)mCalloc(MAX_N_CAL,sizeof(t_cal *));
532
533 Init_Node_Light(n,num);
534
535 return n;
536 }
537
538 //////////////////////////////////////////////////////////////
539 //////////////////////////////////////////////////////////////
540
Make_Node_Lk(t_node * n)541 void Make_Node_Lk(t_node *n)
542 {
543 /* n->n_ex_nodes = (int *)mCalloc(2,sizeof(int)); */
544 return;
545 }
546
547 //////////////////////////////////////////////////////////////
548 //////////////////////////////////////////////////////////////
549
Make_Nexus_Com()550 nexcom **Make_Nexus_Com()
551 {
552 nexcom **com;
553 int i;
554
555 com = (nexcom **)mCalloc(N_MAX_NEX_COM,sizeof(nexcom *));
556
557 for(i=0;i<N_MAX_NEX_COM;i++)
558 {
559 com[i] = (nexcom *)mCalloc(1,sizeof(nexcom));
560 com[i]->name = (char *)mCalloc(T_MAX_NEX_COM,sizeof(char));
561 com[i]->parm = (nexparm **)mCalloc(N_MAX_NEX_PARM,sizeof(nexparm *));
562 }
563
564 return com;
565 }
566
567 //////////////////////////////////////////////////////////////
568 //////////////////////////////////////////////////////////////
569
Make_Nexus_Parm()570 nexparm *Make_Nexus_Parm()
571 {
572 nexparm *parm;
573
574 parm = (nexparm *)mCalloc(1,sizeof(nexparm));
575 parm->name = (char *)mCalloc(T_MAX_TOKEN,sizeof(char ));
576 parm->value = (char *)mCalloc(T_MAX_TOKEN,sizeof(char ));
577
578 return parm;
579 }
580
581 //////////////////////////////////////////////////////////////
582 //////////////////////////////////////////////////////////////
583
Make_Mat(int n_otu)584 matrix *Make_Mat(int n_otu)
585 {
586 matrix *mat;
587 int i;
588
589 mat = (matrix *)mCalloc(1,sizeof(matrix));
590
591 mat->n_otu = n_otu;
592
593 mat->P = (phydbl **)mCalloc(n_otu,sizeof(phydbl *));
594 mat->Q = (phydbl **)mCalloc(n_otu,sizeof(phydbl *));
595 mat->dist = (phydbl **)mCalloc(n_otu,sizeof(phydbl *));
596 mat->on_off = (int *)mCalloc(n_otu,sizeof(int));
597 mat->name = (char **)mCalloc(n_otu,sizeof(char *));
598 mat->tip_node = (t_node **)mCalloc(n_otu,sizeof(t_node *));
599
600
601 for(i=0;i<n_otu;i++)
602 {
603 mat->P[i] = (phydbl *)mCalloc(n_otu,sizeof(phydbl));
604 mat->Q[i] = (phydbl *)mCalloc(n_otu,sizeof(phydbl));
605 mat->dist[i] = (phydbl *)mCalloc(n_otu,sizeof(phydbl));
606 mat->name[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char));
607 }
608
609 return mat;
610 }
611
612 //////////////////////////////////////////////////////////////
613 //////////////////////////////////////////////////////////////
614
Make_Tree_From_Scratch(int n_otu,calign * data)615 t_tree *Make_Tree_From_Scratch(int n_otu, calign *data)
616 {
617 t_tree *tree;
618
619 tree = Make_Tree(n_otu);
620 Init_Tree(tree,n_otu);
621 Make_All_Tree_Nodes(tree);
622 Make_All_Tree_Edges(tree);
623 Make_Tree_Path(tree);
624 if(data)
625 {
626 Copy_Tax_Names_To_Tip_Labels(tree,data);
627 tree->data = data;
628 }
629
630
631 #ifdef BEAGLE
632 //offset the branch's partial indices because BEAGLE insists on first storing the tips/taxa
633 int num_branches = 2*tree->n_otu-1;
634 int i;
635 for(i=0;i<2*tree->n_otu-1;++i)
636 {
637 //For edgeX, its "left" partial lies at index `num_tax + edgeX->num"
638 tree->a_edges[i]->p_lk_left_idx = tree->n_otu + tree->a_edges[i]->p_lk_left_idx;
639 //For edgeX, its "right" partial lies at index `num_tax + edgeX->num + num_branches"
640 tree->a_edges[i]->p_lk_rght_idx = tree->n_otu + tree->a_edges[i]->p_lk_left_idx + num_branches;
641 }
642 #endif
643
644 return tree;
645 }
646
647 //////////////////////////////////////////////////////////////
648 //////////////////////////////////////////////////////////////
649
Make_Tree(int n_otu)650 t_tree *Make_Tree(int n_otu)
651 {
652 t_tree *tree;
653 tree = (t_tree *)mCalloc(1,sizeof(t_tree ));
654 tree->t_dir = (short int *)mCalloc((2*n_otu-2)*(2*n_otu-2),sizeof(int));
655 return tree;
656 }
657
658 //////////////////////////////////////////////////////////////
659 //////////////////////////////////////////////////////////////
660
661
Make_Tree_Path(t_tree * tree)662 void Make_Tree_Path(t_tree *tree)
663 {
664 tree->curr_path = (t_node **)mCalloc(tree->n_otu,sizeof(t_node *));
665 }
666
667 //////////////////////////////////////////////////////////////
668 //////////////////////////////////////////////////////////////
669
670
Make_All_Tree_Nodes(t_tree * tree)671 void Make_All_Tree_Nodes(t_tree *tree)
672 {
673 int i;
674
675 tree->a_nodes = (t_node **)mCalloc(2*tree->n_otu-1,sizeof(t_node *));
676
677 for(i=0;i<2*tree->n_otu-1;++i)
678 {
679 tree->a_nodes[i] = (t_node *)Make_Node_Light(i);
680 if(i < tree->n_otu) tree->a_nodes[i]->tax = YES;
681 else tree->a_nodes[i]->tax = NO;
682 }
683 }
684
685 //////////////////////////////////////////////////////////////
686 //////////////////////////////////////////////////////////////
687
Make_All_Tree_Edges(t_tree * tree)688 void Make_All_Tree_Edges(t_tree *tree)
689 {
690 int i;
691 tree->a_edges = (t_edge **)mCalloc(2*tree->n_otu-1,sizeof(t_edge *));
692 for(i=0;i<2*tree->n_otu-1;++i) tree->a_edges[i] = (t_edge *)Make_Edge_Light(NULL,NULL,i);
693 }
694
695 //////////////////////////////////////////////////////////////
696 //////////////////////////////////////////////////////////////
697
Make_Calign(int n_otu,int crunch_len,int state_len,int init_len,char ** sp_names_in,int n_rm,char ** sp_names_out)698 calign *Make_Calign(int n_otu, int crunch_len, int state_len, int init_len, char **sp_names_in, int n_rm, char **sp_names_out)
699 {
700 calign *cdata;
701 int j;
702
703 cdata = (calign *)mCalloc(1,sizeof(calign));
704 cdata->c_seq = (align **)mCalloc(n_otu,sizeof(align *));
705 cdata->obs_state_frq = (phydbl *)mCalloc(T_MAX_ALPHABET,sizeof(phydbl));
706 cdata->wght = (phydbl *)mCalloc(crunch_len,sizeof(phydbl));
707 cdata->ambigu = (short int *)mCalloc(crunch_len,sizeof(short int));
708 cdata->invar = (short int *)mCalloc(crunch_len,sizeof(short int));
709 cdata->sitepatt = (int *)mCalloc(init_len,sizeof(int ));
710
711 if(n_rm > 0) cdata->c_seq_rm = (align **)mCalloc(n_rm,sizeof(align *));
712
713 for(j=0;j<n_otu;j++)
714 {
715 cdata->c_seq[j] = (align *)mCalloc(1,sizeof(align));
716 cdata->c_seq[j]->name = (char *)mCalloc((int)(strlen(sp_names_in[j])+1),sizeof(char));
717 strcpy(cdata->c_seq[j]->name,sp_names_in[j]);
718 cdata->c_seq[j]->state = (char *)mCalloc(crunch_len*state_len+1,sizeof(char));
719 cdata->c_seq[j]->d_state = (short int *)mCalloc(crunch_len*state_len,sizeof(short int));
720 cdata->c_seq[j]->is_ambigu = (short int *)mCalloc(crunch_len,sizeof(short int));
721 }
722
723 for(j=0;j<n_rm;j++)
724 {
725 cdata->c_seq_rm[j] = (align *)mCalloc(1,sizeof(align));
726 cdata->c_seq_rm[j]->name = (char *)mCalloc((int)(strlen(sp_names_out[j])+1),sizeof(char));
727 strcpy(cdata->c_seq_rm[j]->name,sp_names_out[j]);
728 cdata->c_seq_rm[j]->state = (char *)mCalloc(crunch_len*state_len+1,sizeof(char));
729 cdata->c_seq_rm[j]->d_state = (short int *)mCalloc(crunch_len*state_len,sizeof(short int));
730 cdata->c_seq_rm[j]->is_ambigu = (short int *)mCalloc(crunch_len,sizeof(short int));
731 }
732
733 return cdata;
734 }
735
736 //////////////////////////////////////////////////////////////
737 //////////////////////////////////////////////////////////////
738
739
Make_Treelist(int list_size)740 t_treelist *Make_Treelist(int list_size)
741 {
742 t_treelist *tlist;
743
744 tlist = (t_treelist *)mCalloc(1,sizeof(t_treelist));
745 tlist->list_size = list_size;
746 tlist->tree = (t_tree **)mCalloc(list_size,sizeof(t_tree *));
747
748 return tlist;
749 }
750
751 //////////////////////////////////////////////////////////////
752 //////////////////////////////////////////////////////////////
753
Make_Optimiz()754 t_opt *Make_Optimiz()
755 {
756 t_opt *s_opt;
757 s_opt = (t_opt *)mCalloc(1,sizeof(t_opt));
758 return s_opt;
759 }
760
761 //////////////////////////////////////////////////////////////
762 //////////////////////////////////////////////////////////////
763
Make_Custom_Model(t_mod * mod)764 void Make_Custom_Model(t_mod *mod)
765 {
766 if(!mod->r_mat)
767 {
768 PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
769 Warn_And_Exit("");
770 }
771
772 if(!mod->r_mat->rr->v)
773 mod->r_mat->rr->v = (phydbl *)mCalloc(mod->ns*(mod->ns-1)/2,sizeof(phydbl));
774
775 if(!mod->r_mat->rr_val->v)
776 mod->r_mat->rr_val->v = (phydbl *)mCalloc(mod->ns*(mod->ns-1)/2,sizeof(phydbl));
777
778 if(!mod->r_mat->rr_num->v)
779 mod->r_mat->rr_num->v = (int *)mCalloc(mod->ns*(mod->ns-1)/2,sizeof(int *));
780
781 if(!mod->r_mat->n_rr_per_cat->v)
782 mod->r_mat->n_rr_per_cat->v = (int *)mCalloc(mod->ns*(mod->ns-1)/2,sizeof(int));
783 }
784
785 //////////////////////////////////////////////////////////////
786 //////////////////////////////////////////////////////////////
787
Make_String(int len)788 t_string *Make_String(int len)
789 {
790 t_string *ts;
791
792 ts = (t_string *)mCalloc(1,sizeof(t_string));
793 ts->s = (char *)mCalloc(len,sizeof(char));
794
795 return(ts);
796 }
797
798 //////////////////////////////////////////////////////////////
799 //////////////////////////////////////////////////////////////
800
Make_Model_Basic()801 t_mod *Make_Model_Basic()
802 {
803 t_mod *mod;
804
805 mod = (t_mod *)mCalloc(1,sizeof(t_mod));
806
807 mod->modelname = Make_String(T_MAX_NAME);
808 Init_String(mod->modelname);
809
810 mod->custom_mod_string = Make_String(T_MAX_NAME);
811 Init_String(mod->custom_mod_string);
812
813 mod->ras = Make_RAS_Basic();
814
815 mod->kappa = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
816 Init_Scalar_Dbl(mod->kappa);
817
818 mod->lambda = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
819 Init_Scalar_Dbl(mod->lambda);
820
821 mod->br_len_mult = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
822 Init_Scalar_Dbl(mod->br_len_mult);
823
824 mod->br_len_mult_unscaled = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
825 Init_Scalar_Dbl(mod->br_len_mult_unscaled);
826
827 mod->mr = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
828 Init_Scalar_Dbl(mod->mr);
829
830 mod->e_frq_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
831 Init_Scalar_Dbl(mod->e_frq_weight);
832
833 mod->r_mat_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
834 Init_Scalar_Dbl(mod->r_mat_weight);
835
836 mod->aa_rate_mat_file = Make_String(T_MAX_FILE);
837 Init_String(mod->aa_rate_mat_file);
838
839 return mod;
840 }
841
842 //////////////////////////////////////////////////////////////
843 //////////////////////////////////////////////////////////////
844
845 /*! Call only when the values of mod->ns & ras->n_catg is set to its final value */
846
Make_Model_Complete(t_mod * mod)847 void Make_Model_Complete(t_mod *mod)
848 {
849 mod->Pij_rr = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
850 Init_Vect_Dbl(0,mod->Pij_rr);
851
852 #if (defined(__AVX__) || defined(__SSE3__))
853 #ifndef WIN32
854 if(posix_memalign((void **)&mod->Pij_rr->v,BYTE_ALIGN,(size_t)mod->ras->n_catg*mod->ns*mod->ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
855 #else
856 mod->Pij_rr->v = _aligned_malloc(mod->ras->n_catg*mod->ns*mod->ns*sizeof(phydbl),BYTE_ALIGN);
857 #endif
858 #else
859 mod->Pij_rr->v = (phydbl *)mCalloc(mod->ras->n_catg*mod->ns*mod->ns,sizeof(phydbl));
860 #endif
861
862 mod->eigen = (eigen *)Make_Eigen_Struct(mod->ns);
863
864 // If r_mat (e_frq) are not NULL, then they have been created elsewhere and affected.
865 if(!mod->r_mat)
866 {
867 mod->r_mat = (t_rmat *)Make_Rmat(mod->ns);
868 Init_Rmat(mod->r_mat);
869 }
870
871 if(!mod->e_frq)
872 {
873 mod->e_frq = (t_efrq *)Make_Efrq(mod->ns);
874 Init_Efrq(NULL,mod->e_frq);
875 }
876
877 Make_RAS_Complete(mod->ras);
878
879 mod->e_frq->user_b_freq->len = mod->ns;
880
881 if(mod->whichmodel < 0)
882 {
883 PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
884 Exit("\n");
885 }
886
887 if(mod->whichmodel == CUSTOM)
888 {
889 Make_Custom_Model(mod);
890 Translate_Custom_Mod_String(mod);
891 }
892
893 if(mod->io->datatype == NT && mod->whichmodel == GTR)
894 {
895 Make_Custom_Model(mod);
896 }
897 else if(mod->io->datatype == GENERIC)
898 {
899 Make_Custom_Model(mod);
900 }
901 }
902
903 //////////////////////////////////////////////////////////////
904 //////////////////////////////////////////////////////////////
905
Make_RAS_Basic()906 t_ras *Make_RAS_Basic()
907 {
908 t_ras *ras;
909
910 ras = (t_ras *)mCalloc(1,sizeof(t_ras));
911
912 ras->gamma_r_proba = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
913 Init_Vect_Dbl(0,ras->gamma_r_proba);
914 ras->gamma_r_proba->v = NULL;
915
916 ras->gamma_r_proba_unscaled = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
917 Init_Vect_Dbl(0,ras->gamma_r_proba_unscaled);
918 ras->gamma_r_proba_unscaled->v = NULL;
919
920 ras->gamma_rr = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
921 Init_Vect_Dbl(0,ras->gamma_rr);
922 ras->gamma_rr->v = NULL;
923
924 ras->gamma_rr_unscaled = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
925 Init_Vect_Dbl(0,ras->gamma_rr_unscaled);
926 ras->gamma_rr_unscaled->v = NULL;
927
928 ras->alpha = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
929 Init_Scalar_Dbl(ras->alpha);
930
931 ras->pinvar = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
932 Init_Scalar_Dbl(ras->pinvar);
933
934 ras->free_rate_mr = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
935 Init_Scalar_Dbl(ras->free_rate_mr);
936
937 return(ras);
938 }
939
940 //////////////////////////////////////////////////////////////
941 //////////////////////////////////////////////////////////////
942
943 /*! Call only when the value of ras->n_catg is set to its final value */
Make_RAS_Complete(t_ras * ras)944 void Make_RAS_Complete(t_ras *ras)
945 {
946 if(!ras->gamma_r_proba->v)
947 {
948 ras->gamma_r_proba->v = (phydbl *)mCalloc(ras->n_catg,sizeof(phydbl));
949 ras->gamma_r_proba_unscaled->v = (phydbl *)mCalloc(ras->n_catg,sizeof(phydbl));
950 ras->gamma_rr->v = (phydbl *)mCalloc(ras->n_catg,sizeof(phydbl));
951 ras->gamma_rr_unscaled->v = (phydbl *)mCalloc(ras->n_catg,sizeof(phydbl));
952 ras->skip_rate_cat = (short int *)mCalloc(ras->n_catg,sizeof(short int));
953 }
954 }
955
956 //////////////////////////////////////////////////////////////
957 //////////////////////////////////////////////////////////////
958
Make_Efrq(int ns)959 t_efrq *Make_Efrq(int ns)
960 {
961 t_efrq *e_frq;
962
963 e_frq = (t_efrq *)mCalloc(1,sizeof(t_efrq));
964
965 e_frq->pi = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
966 e_frq->pi->len = ns;
967
968 #if (defined(__AVX__) || defined(__SSE3__))
969 #ifndef WIN32
970 if(posix_memalign((void **)&e_frq->pi->v,BYTE_ALIGN,(size_t)ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
971 #else
972 e_frq->pi->v = _aligned_malloc(ns * sizeof(phydbl),BYTE_ALIGN);
973 #endif
974 #else
975 e_frq->pi->v = (phydbl *)mCalloc(ns,sizeof(phydbl));
976 #endif
977
978
979 e_frq->pi_unscaled = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
980 e_frq->pi_unscaled->v = (phydbl *)mCalloc(ns,sizeof(phydbl));
981 e_frq->pi_unscaled->len = ns;
982
983 e_frq->user_b_freq = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
984 Init_Vect_Dbl(0,e_frq->user_b_freq);
985 e_frq->user_b_freq->v = (phydbl *)mCalloc(T_MAX_OPTION,sizeof(phydbl));
986
987 return e_frq;
988 }
989
990 //////////////////////////////////////////////////////////////
991 //////////////////////////////////////////////////////////////
992
Make_Rmat(int ns)993 t_rmat *Make_Rmat(int ns)
994 {
995 t_rmat *r_mat;
996
997 r_mat = (t_rmat *)mCalloc(1,sizeof(t_rmat));
998
999 r_mat->qmat = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
1000 Init_Vect_Dbl(0,r_mat->qmat);
1001
1002 r_mat->qmat_buff = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
1003 Init_Vect_Dbl(0,r_mat->qmat_buff);
1004
1005 r_mat->rr = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
1006 Init_Vect_Dbl(0,r_mat->rr);
1007
1008 r_mat->rr_val = (vect_dbl *)mCalloc(1,sizeof(vect_dbl));
1009 Init_Vect_Dbl(0,r_mat->rr_val);
1010
1011 r_mat->rr_num = (vect_int *)mCalloc(1,sizeof(vect_int));
1012 Init_Vect_Int(0,r_mat->rr_num);
1013
1014 r_mat->n_rr_per_cat = (vect_int *)mCalloc(1,sizeof(vect_int));
1015 Init_Vect_Int(0,r_mat->n_rr_per_cat);
1016
1017 r_mat->qmat->v = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
1018 r_mat->qmat_buff->v = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
1019
1020 return(r_mat);
1021 }
1022
1023 //////////////////////////////////////////////////////////////
1024 //////////////////////////////////////////////////////////////
1025
Make_Input()1026 option *Make_Input()
1027 {
1028 int i;
1029 option* io = (option *)mCalloc(1,sizeof(option));
1030
1031 io->in_align_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1032 io->in_tree_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1033 io->in_constraint_tree_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1034 io->in_coord_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1035 io->out_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1036 io->out_tree_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1037 io->out_trees_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1038 io->out_ancestral_seq_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1039 io->out_ancestral_tree_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1040 io->out_boot_tree_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1041 io->out_boot_stats_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1042 io->out_stats_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1043 io->out_lk_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1044 io->weight_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1045 io->out_ps_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1046 io->out_summary_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1047 io->out_trace_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1048 io->out_json_trace_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1049 io->nt_or_cd = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1050 io->run_id_string = (char *)mCalloc(T_MAX_OPTION,sizeof(char));
1051 io->clade_list_file = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1052 io->alphabet = (char **)mCalloc(T_MAX_ALPHABET,sizeof(char *));
1053 for(i=0;i<T_MAX_ALPHABET;i++) io->alphabet[i] = (char *)mCalloc(T_MAX_STATE,sizeof(char ));
1054 io->treelist = (t_treelist *)mCalloc(1,sizeof(t_treelist));
1055 io->mcmc = (t_mcmc *)MCMC_Make_MCMC_Struct();
1056 io->rates = (t_rate *)RATES_Make_Rate_Struct(-1);
1057 io->times = (t_time *)TIMES_Make_Time_Struct(-1);
1058
1059 return io;
1060 }
1061
1062 //////////////////////////////////////////////////////////////
1063 //////////////////////////////////////////////////////////////
1064
MCMC_Make_MCMC_Struct()1065 t_mcmc *MCMC_Make_MCMC_Struct()
1066 {
1067 t_mcmc *mcmc;
1068
1069 mcmc = (t_mcmc *)mCalloc(1,sizeof(t_mcmc));
1070 mcmc->out_filename = (char *)mCalloc(T_MAX_FILE,sizeof(char));
1071
1072 return(mcmc);
1073 }
1074
1075 //////////////////////////////////////////////////////////////
1076 //////////////////////////////////////////////////////////////
1077
Make_Eigen_Struct(int ns)1078 eigen *Make_Eigen_Struct(int ns)
1079 {
1080 eigen *eig;
1081
1082 eig = (eigen *)mCalloc(1,sizeof(eigen));
1083 eig->size = ns;
1084 eig->space = (phydbl *)mCalloc(2*ns,sizeof(phydbl));
1085 eig->space_int = (int *)mCalloc(2*ns,sizeof(int));
1086 eig->e_val = (phydbl *)mCalloc(ns,sizeof(phydbl));
1087 eig->e_val_im = (phydbl *)mCalloc(ns,sizeof(phydbl));
1088 eig->r_e_vect_im = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
1089 eig->q = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
1090
1091 #if (defined(__AVX__) || defined(__SSE3__))
1092 #ifndef WIN32
1093 if(posix_memalign((void **)&eig->r_e_vect,BYTE_ALIGN,ns*ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
1094 if(posix_memalign((void **)&eig->l_e_vect,BYTE_ALIGN,ns*ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
1095 if(posix_memalign((void **)&eig->dum,BYTE_ALIGN,ns*ns*sizeof(phydbl))) Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
1096 #else
1097 eig->r_e_vect = _aligned_malloc(ns*ns*sizeof(phydbl),BYTE_ALIGN);
1098 eig->l_e_vect = _aligned_malloc(ns*ns*sizeof(phydbl),BYTE_ALIGN);
1099 eig->dum = _aligned_malloc(ns*ns*sizeof(phydbl),BYTE_ALIGN);
1100 #endif
1101 #else
1102 eig->r_e_vect = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
1103 eig->l_e_vect = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
1104 eig->dum = (phydbl *)mCalloc(ns*ns,sizeof(phydbl));
1105 #endif
1106
1107 Init_Eigen_Struct(eig);
1108
1109 return eig;
1110 }
1111
1112 //////////////////////////////////////////////////////////////
1113 //////////////////////////////////////////////////////////////
1114
1115 //////////////////////////////////////////////////////////////
1116 //////////////////////////////////////////////////////////////
1117
Make_Short_L(t_tree * tree)1118 void Make_Short_L(t_tree *tree)
1119 {
1120 if(!tree->short_l)
1121 tree->short_l = (phydbl *)mCalloc(tree->n_short_l,sizeof(phydbl));
1122 }
1123
1124 //////////////////////////////////////////////////////////////
1125 //////////////////////////////////////////////////////////////
1126
XML_Make_Attribute(xml_attr * prev,char * attr_name,char * attr_value)1127 xml_attr *XML_Make_Attribute(xml_attr *prev, char *attr_name, char *attr_value)
1128 {
1129 xml_attr *new_attr;
1130
1131 new_attr = (xml_attr *)mCalloc(1,sizeof(xml_attr));
1132
1133 new_attr->prev = prev;
1134 new_attr->next = NULL;
1135 if(prev != NULL) prev->next = new_attr;
1136
1137 new_attr->name = (char *)mCalloc(strlen(attr_name)+1,sizeof(char));
1138 strcpy(new_attr->name,attr_name);
1139
1140 new_attr->value = (char *)mCalloc(strlen(attr_value)+1,sizeof(char));
1141 strcpy(new_attr->value,attr_value);
1142
1143 return new_attr;
1144 }
1145
1146 //////////////////////////////////////////////////////////////
1147 //////////////////////////////////////////////////////////////
1148
XML_Make_Node_Id(xml_node * n,char * id)1149 void XML_Make_Node_Id(xml_node *n, char *id)
1150 {
1151 if(id) n->id = (char *)mCalloc(strlen(id)+1,sizeof(char));
1152 }
1153
1154 //////////////////////////////////////////////////////////////
1155 //////////////////////////////////////////////////////////////
1156
XML_Make_Node_Value(xml_node * n,char * val)1157 void XML_Make_Node_Value(xml_node *n, char *val)
1158 {
1159 if(val) n->value = (char *)mCalloc(strlen(val)+1,sizeof(char));
1160 }
1161
1162 //////////////////////////////////////////////////////////////
1163 //////////////////////////////////////////////////////////////
1164
1165
XML_Make_Node(char * name)1166 xml_node *XML_Make_Node(char *name)
1167 {
1168 xml_node *new_node = (xml_node *)mCalloc(1,sizeof(xml_node));
1169 if(name) new_node->name = (char *)mCalloc(strlen(name)+1,sizeof(char));
1170 new_node->ds = (t_ds *)mCalloc(1,sizeof(t_ds));
1171 return new_node;
1172 }
1173
1174 //////////////////////////////////////////////////////////////
1175 //////////////////////////////////////////////////////////////
1176
Make_Best_Spr(t_tree * tree)1177 void Make_Best_Spr(t_tree *tree)
1178 {
1179 tree->best_spr = Make_One_Spr(tree);
1180 Init_One_Spr(tree->best_spr);
1181 }
1182
1183 //////////////////////////////////////////////////////////////
1184 //////////////////////////////////////////////////////////////
1185
Make_Spr(t_tree * tree)1186 void Make_Spr(t_tree *tree)
1187 {
1188 Make_Spr_List_One_Edge(tree);
1189 Make_Spr_List_All_Edge(tree);
1190 Make_Best_Spr(tree);
1191 }
1192
1193 //////////////////////////////////////////////////////////////
1194 //////////////////////////////////////////////////////////////
1195
Make_Spr_List_One_Edge(t_tree * tree)1196 void Make_Spr_List_One_Edge(t_tree *tree)
1197 {
1198 int i;
1199
1200 tree->size_spr_list_one_edge = 2*tree->n_otu-3;
1201
1202 tree->spr_list_one_edge = (t_spr **)mCalloc(2*tree->n_otu-2,sizeof(t_spr *));
1203
1204 for(i=0;i<2*tree->n_otu-2;++i)
1205 {
1206 tree->spr_list_one_edge[i] = Make_One_Spr(tree);
1207 Init_One_Spr(tree->spr_list_one_edge[i]);
1208 }
1209
1210 tree->perform_spr_right_away = NO;
1211 }
1212
1213 //////////////////////////////////////////////////////////////
1214 //////////////////////////////////////////////////////////////
1215
Make_Spr_List_All_Edge(t_tree * tree)1216 void Make_Spr_List_All_Edge(t_tree *tree)
1217 {
1218 int i;
1219
1220 tree->size_spr_list_all_edge = 2*tree->n_otu-3;
1221
1222 tree->spr_list_all_edge = (t_spr **)mCalloc(2*tree->n_otu-2,sizeof(t_spr *));
1223
1224 for(i=0;i<2*tree->n_otu-2;++i)
1225 {
1226 tree->spr_list_all_edge[i] = Make_One_Spr(tree);
1227 Init_One_Spr(tree->spr_list_all_edge[i]);
1228 }
1229
1230 tree->perform_spr_right_away = NO;
1231 }
1232
1233 //////////////////////////////////////////////////////////////
1234 //////////////////////////////////////////////////////////////
1235
Make_One_Spr(t_tree * tree)1236 t_spr *Make_One_Spr(t_tree *tree)
1237 {
1238 t_spr *a_spr;
1239 a_spr = (t_spr *)mCalloc(1,sizeof(t_spr));
1240 a_spr->path = (t_node **)mCalloc(tree->n_otu,sizeof(t_node *));
1241 return a_spr;
1242 }
1243
1244 //////////////////////////////////////////////////////////////
1245 //////////////////////////////////////////////////////////////
1246
TIMES_Make_Time_Struct(int n_otu)1247 t_time *TIMES_Make_Time_Struct(int n_otu)
1248 {
1249 t_time *times;
1250
1251 times = (t_time *)mCalloc(1,sizeof(t_time));
1252
1253 if(n_otu > 0)
1254 {
1255 times->calib_prob = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1256 times->t_prior_min_ori = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1257 times->t_prior_max_ori = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1258 times->times_partial_proba = (phydbl *)mCalloc(n_otu*n_otu,sizeof(phydbl));
1259 times->numb_calib_chosen = (int *)mCalloc(n_otu*n_otu,sizeof(phydbl));
1260 times->a_cal = (t_cal **)mCalloc(10*n_otu,sizeof(t_cal *));
1261
1262 times->nd_t = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1263 times->buff_t = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1264 times->true_t = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1265 times->t_mean = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1266 times->t_prior = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1267 times->t_prior_min = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1268 times->t_prior_max = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1269 times->t_floor = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1270 times->t_rank = (int *)mCalloc(2*n_otu-1,sizeof(int));
1271 times->t_has_prior = (short int *)mCalloc(2*n_otu-1,sizeof(short int));
1272
1273 times->mean_t = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1274
1275 times->n_jps = (int *)mCalloc(2*n_otu-1,sizeof(int));
1276 times->t_jps = (int *)mCalloc(2*n_otu-2,sizeof(int));
1277
1278 times->time_slice_lims = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1279 times->n_time_slice_spans = (int *)mCalloc(2*n_otu-1,sizeof(int));
1280 times->curr_slice = (int *)mCalloc(2*n_otu-1,sizeof(int));
1281 times->has_survived = (int *)mCalloc(2*n_otu-1,sizeof(int));
1282 }
1283 return(times);
1284 }
1285
1286 //////////////////////////////////////////////////////////////
1287 //////////////////////////////////////////////////////////////
1288
RATES_Make_Rate_Struct(int n_otu)1289 t_rate *RATES_Make_Rate_Struct(int n_otu)
1290 {
1291 t_rate *rates;
1292
1293 rates = (t_rate *)mCalloc(1,sizeof(t_rate));
1294 rates->is_allocated = NO;
1295
1296 if(n_otu > 0)
1297 {
1298 rates->is_allocated = YES;
1299 rates->nd_r = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1300 rates->br_r = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1301 rates->buff_br_r = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1302 rates->buff_nd_r = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1303 rates->true_r = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1304 rates->dens = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl));
1305 rates->triplet = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1306 rates->cov_l = (phydbl *)mCalloc((2*n_otu-2)*(2*n_otu-2),sizeof(phydbl));
1307 rates->invcov = (phydbl *)mCalloc((2*n_otu-2)*(2*n_otu-2),sizeof(phydbl));
1308 rates->mean_l = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl));
1309 rates->ml_l = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl));
1310 rates->cur_l = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl));
1311 rates->u_ml_l = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1312 rates->u_cur_l = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1313 rates->cov_r = (phydbl *)mCalloc((2*n_otu-2)*(2*n_otu-2),sizeof(phydbl));
1314 rates->cond_var = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl));
1315 rates->mean_r = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl));
1316 rates->lca = (t_node **)mCalloc((2*n_otu-1)*(2*n_otu-1),sizeof(t_node *));
1317 rates->reg_coeff = (phydbl *)mCalloc((2*n_otu-2)*(2*n_otu-2),sizeof(phydbl));
1318 rates->trip_reg_coeff = (phydbl *)mCalloc((2*n_otu-2)*(6*n_otu-9),sizeof(phydbl));
1319 rates->trip_cond_cov = (phydbl *)mCalloc((2*n_otu-2)*9,sizeof(phydbl));
1320 rates->_2n_vect1 = (phydbl *)mCalloc(2*n_otu,sizeof(phydbl));
1321 rates->_2n_vect2 = (phydbl *)mCalloc(2*n_otu,sizeof(phydbl));
1322 rates->_2n_vect3 = (phydbl *)mCalloc(2*n_otu,sizeof(phydbl));
1323 rates->_2n_vect4 = (phydbl *)mCalloc(2*n_otu,sizeof(phydbl));
1324 rates->_2n_vect5 = (short int *)mCalloc(2*n_otu,sizeof(short int));
1325 rates->_2n2n_vect1 = (phydbl *)mCalloc(4*n_otu*n_otu,sizeof(phydbl));
1326 rates->_2n2n_vect2 = (phydbl *)mCalloc(4*n_otu*n_otu,sizeof(phydbl));
1327 rates->br_do_updt = (short int *)mCalloc(2*n_otu-1,sizeof(short int));
1328 rates->cur_gamma_prior_mean = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1329 rates->cur_gamma_prior_var = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl));
1330 rates->n_tips_below = (int *)mCalloc(2*n_otu-1,sizeof(int));
1331 rates->model_name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1332 }
1333
1334 return rates;
1335 }
1336
1337 //////////////////////////////////////////////////////////////
1338 //////////////////////////////////////////////////////////////
1339
Make_Calibration()1340 t_cal *Make_Calibration()
1341 {
1342 t_cal *calib;
1343 calib = (t_cal *)mCalloc(1, sizeof(t_cal));
1344 return(calib);
1345 }
1346
1347 //////////////////////////////////////////////////////////////
1348 //////////////////////////////////////////////////////////////
1349
Make_Clade()1350 t_clad *Make_Clade()
1351 {
1352 t_clad *clade;
1353 clade = (t_clad *)mCalloc(1, sizeof(t_clad));
1354 return(clade);
1355 }
1356
1357 //////////////////////////////////////////////////////////////
1358 //////////////////////////////////////////////////////////////
1359
1360
Make_Target_Tip(int n)1361 t_node **Make_Target_Tip(int n)
1362 {
1363 t_node **this;
1364 this = (t_node **)mCalloc(n,sizeof(t_node *));
1365 return(this);
1366 }
1367
1368 //////////////////////////////////////////////////////////////
1369 //////////////////////////////////////////////////////////////
1370
Make_All_Calibration(t_tree * tree)1371 void Make_All_Calibration(t_tree *tree)
1372 {
1373 int i;
1374 t_cal **all_cal;
1375
1376 assert(tree->rates);
1377
1378 all_cal = (t_cal **)mCalloc(2*tree->n_otu-1,sizeof(t_cal *));
1379
1380 For(i,2*tree->n_otu-1) all_cal[i] = Make_Calibration();
1381
1382 tree->times->a_cal = all_cal;
1383
1384 }
1385
1386 //////////////////////////////////////////////////////////////
1387 //////////////////////////////////////////////////////////////
1388
Make_Rmat_Weight(t_tree * mixt_tree)1389 void Make_Rmat_Weight(t_tree *mixt_tree)
1390 {
1391 t_tree *tree, *buff_tree;
1392 scalar_dbl *curr_weight;
1393
1394 tree = mixt_tree;
1395 do
1396 {
1397 if(tree->is_mixt_tree == YES) tree = tree->next;
1398
1399 buff_tree = mixt_tree->next;
1400 do
1401 {
1402 if(buff_tree->mod->r_mat_weight == tree->mod->r_mat_weight) break;
1403 buff_tree = buff_tree->next;
1404 }
1405 while(buff_tree != tree);
1406
1407 if(buff_tree == tree) Free(tree->mod->r_mat_weight);
1408
1409 tree = tree->next;
1410 }
1411 while(tree);
1412
1413
1414 tree = mixt_tree;
1415 do
1416 {
1417 if(tree->is_mixt_tree == YES) tree = tree->next;
1418 tree->mod->r_mat_weight = NULL;
1419 tree = tree->next;
1420 }
1421 while(tree);
1422
1423
1424 tree = mixt_tree->next;
1425 tree->mod->r_mat_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
1426 Init_Scalar_Dbl(tree->mod->r_mat_weight);
1427 tree->mod->r_mat_weight->v = 1.0;
1428 curr_weight = tree->mod->r_mat_weight;
1429
1430
1431 buff_tree = tree = mixt_tree;
1432 do // For each mixt_tree
1433 {
1434 if(tree->is_mixt_tree == YES)
1435 {
1436 tree = tree->next;
1437 }
1438
1439 buff_tree = mixt_tree->next;
1440 do
1441 {
1442 if(buff_tree->mod->r_mat == tree->mod->r_mat)
1443 {
1444 tree->mod->r_mat_weight = buff_tree->mod->r_mat_weight;
1445 break;
1446 }
1447 buff_tree = buff_tree->next;
1448 }
1449 while(buff_tree != tree);
1450
1451 if(!tree->mod->r_mat_weight)
1452 {
1453 tree->mod->r_mat_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
1454 Init_Scalar_Dbl(tree->mod->r_mat_weight);
1455 tree->mod->r_mat_weight->v = 1.0;
1456 curr_weight->next = tree->mod->r_mat_weight;
1457 tree->mod->r_mat_weight->prev = curr_weight;
1458 curr_weight = tree->mod->r_mat_weight;
1459 }
1460
1461 tree = tree->next;
1462
1463 }
1464 while(tree);
1465
1466 }
1467
1468 //////////////////////////////////////////////////////////////
1469 //////////////////////////////////////////////////////////////
1470
Make_Efrq_Weight(t_tree * mixt_tree)1471 void Make_Efrq_Weight(t_tree *mixt_tree)
1472 {
1473 t_tree *tree, *buff_tree;
1474 scalar_dbl *curr_weight;
1475
1476
1477 tree = mixt_tree;
1478 do
1479 {
1480 if(tree->is_mixt_tree == YES) tree = tree->next;
1481
1482 buff_tree = mixt_tree->next;
1483 do
1484 {
1485 if(buff_tree->mod->e_frq_weight == tree->mod->e_frq_weight) break;
1486 buff_tree = buff_tree->next;
1487 }
1488 while(buff_tree != tree);
1489
1490 if(buff_tree == tree)
1491 {
1492 Free(tree->mod->e_frq_weight);
1493 }
1494
1495 tree = tree->next;
1496 }
1497 while(tree);
1498
1499
1500 tree = mixt_tree;
1501 do
1502 {
1503 if(tree->is_mixt_tree == YES) tree = tree->next;
1504 tree->mod->e_frq_weight = NULL;
1505 tree = tree->next;
1506 }
1507 while(tree);
1508
1509
1510 tree = mixt_tree->next;
1511 tree->mod->e_frq_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
1512 Init_Scalar_Dbl(tree->mod->e_frq_weight);
1513 tree->mod->e_frq_weight->v = 1.0;
1514 curr_weight = tree->mod->e_frq_weight;
1515
1516
1517 buff_tree = tree = mixt_tree;
1518 do // For each mixt_tree
1519 {
1520 if(tree->is_mixt_tree == YES) tree = tree->next;
1521
1522 buff_tree = mixt_tree->next;
1523 do
1524 {
1525 if(buff_tree->mod->e_frq == tree->mod->e_frq)
1526 {
1527 tree->mod->e_frq_weight = buff_tree->mod->e_frq_weight;
1528 break;
1529 }
1530 buff_tree = buff_tree->next;
1531 }
1532 while(buff_tree != tree);
1533
1534 if(!tree->mod->e_frq_weight)
1535 {
1536 tree->mod->e_frq_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl));
1537 Init_Scalar_Dbl(tree->mod->e_frq_weight);
1538 tree->mod->e_frq_weight->v = 1.0;
1539 curr_weight->next = tree->mod->e_frq_weight;
1540 tree->mod->e_frq_weight->prev = curr_weight;
1541 curr_weight = tree->mod->e_frq_weight;
1542 }
1543
1544 tree = tree->next;
1545
1546 }
1547 while(tree);
1548 }
1549
1550 //////////////////////////////////////////////////////////////
1551 //////////////////////////////////////////////////////////////
1552
1553 //////////////////////////////////////////////////////////////
1554 //////////////////////////////////////////////////////////////
1555
GEO_Make_Geo_Basic()1556 t_geo *GEO_Make_Geo_Basic()
1557 {
1558 t_geo *t;
1559 t = (t_geo *)mCalloc(1,sizeof(t_geo));
1560 return(t);
1561 }
1562
1563 //////////////////////////////////////////////////////////////
1564 //////////////////////////////////////////////////////////////
1565
GEO_Make_Geo_Complete(int ldscape_sz,int n_dim,int n_tax,t_geo * t)1566 void GEO_Make_Geo_Complete(int ldscape_sz, int n_dim, int n_tax, t_geo *t)
1567 {
1568 int i;
1569
1570 // F matrix
1571 t->f_mat = (phydbl *)mCalloc(ldscape_sz*ldscape_sz,sizeof(phydbl));
1572
1573 // R matrix
1574 t->r_mat = (phydbl *)mCalloc(ldscape_sz*ldscape_sz,sizeof(phydbl));
1575
1576 // Occupation vectors: one vector for each node
1577 t->occup = (int *)mCalloc((2*n_tax-1)*ldscape_sz,sizeof(int));
1578
1579 // Lineage locations
1580 t->idx_loc = (int *)mCalloc((int)(2*n_tax-1),sizeof(int));
1581
1582 // Sorted node heights
1583 t->sorted_nd = (t_node **)mCalloc((int)(2*n_tax-1),sizeof(t_node *));
1584
1585 // Covariance matrix
1586 t->cov = (phydbl *)mCalloc((int)(n_dim*n_dim),sizeof(phydbl));
1587
1588 // gives the location occupied beneath each node in the tree
1589 t->idx_loc_beneath = (int *)mCalloc((int)(2*n_tax-1)*ldscape_sz,sizeof(int));
1590
1591 // Locations
1592 t->coord_loc = (t_geo_coord **)mCalloc(ldscape_sz,sizeof(t_geo_coord *));
1593 for(i=0;i<ldscape_sz;i++) t->coord_loc[i] = GEO_Make_Geo_Coord(n_dim);
1594 }
1595
1596 //////////////////////////////////////////////////////////////
1597 //////////////////////////////////////////////////////////////
1598
GEO_Make_Geo_Coord(int dim)1599 t_geo_coord *GEO_Make_Geo_Coord(int dim)
1600 {
1601 t_geo_coord *t;
1602 t = (t_geo_coord *)mCalloc(1,sizeof(t_geo_coord));
1603 t->lonlat = (phydbl *)mCalloc(dim,sizeof(phydbl));
1604 t->id = (char *)mCalloc(T_MAX_ID_COORD,sizeof(char));
1605
1606 t->cpy = (t_geo_coord *)mCalloc(1,sizeof(t_geo_coord));
1607 t->cpy->lonlat = (phydbl *)mCalloc(dim,sizeof(phydbl));
1608 t->cpy->id = (char *)mCalloc(T_MAX_ID_COORD,sizeof(char));
1609
1610 return(t);
1611 }
1612
1613 //////////////////////////////////////////////////////////////
1614 //////////////////////////////////////////////////////////////
1615
PHYREX_Make_Migrep_Model(int n_otu,int dim)1616 t_phyrex_mod *PHYREX_Make_Migrep_Model(int n_otu,int dim)
1617 {
1618 t_phyrex_mod *t;
1619 t = (t_phyrex_mod *)mCalloc(1,sizeof(t_phyrex_mod));
1620 t->lim_up = GEO_Make_Geo_Coord(dim);
1621 t->lim_do = GEO_Make_Geo_Coord(dim);
1622 t->sigsq_scale = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl));
1623 return(t);
1624 }
1625
1626 //////////////////////////////////////////////////////////////
1627 //////////////////////////////////////////////////////////////
1628
PHYREX_Make_Disk_Event(int n_dim,int n_otu)1629 t_dsk *PHYREX_Make_Disk_Event(int n_dim, int n_otu)
1630 {
1631 t_dsk *t;
1632
1633 t = (t_dsk *)mCalloc(1,sizeof(t_dsk));
1634 t->centr = GEO_Make_Geo_Coord(n_dim);
1635 t->id = (char *)mCalloc(T_MAX_ID_DISK,sizeof(char));
1636 t->ldsk_a = (t_ldsk **)mCalloc(n_otu,sizeof(t_ldsk *));
1637
1638 return(t);
1639 }
1640
1641 //////////////////////////////////////////////////////////////
1642 //////////////////////////////////////////////////////////////
1643
PHYREX_Make_Lindisk_Node(int n_dim)1644 t_ldsk *PHYREX_Make_Lindisk_Node(int n_dim)
1645 {
1646 t_ldsk *t;
1647 t = (t_ldsk *)mCalloc(1,sizeof(t_ldsk));
1648 t->coord = GEO_Make_Geo_Coord(n_dim);
1649 t->cpy_coord = GEO_Make_Geo_Coord(n_dim);
1650 return(t);
1651 }
1652 //////////////////////////////////////////////////////////////
1653 //////////////////////////////////////////////////////////////
1654
PHYREX_Make_Lindisk_Next(t_ldsk * t)1655 void PHYREX_Make_Lindisk_Next(t_ldsk *t)
1656 {
1657
1658 if(t->n_next == 0)
1659 t->next = (t_ldsk **)mCalloc(NEXT_BLOCK_SIZE,sizeof(t_ldsk *));
1660 else if(!(t->n_next%NEXT_BLOCK_SIZE))
1661 t->next = (t_ldsk **)mRealloc(t->next,t->n_next+NEXT_BLOCK_SIZE,sizeof(t_ldsk *));
1662
1663 t->n_next++;
1664 /* PhyML_Printf("\n. make next for ldsk %s n_next set to %d [%p] %d",t->coord->id,t->n_next,t->next,NEXT_BLOCK_SIZE); */
1665 }
1666
1667 //////////////////////////////////////////////////////////////
1668 //////////////////////////////////////////////////////////////
1669
Make_Poly(int n)1670 t_poly *Make_Poly(int n)
1671 {
1672 t_poly *p;
1673 int i;
1674 p = (t_poly *)mCalloc(1,sizeof(t_poly));
1675 p->poly_vert = (t_geo_coord **)mCalloc(n,sizeof(t_geo_coord *));
1676 for(i=0;i<n;i++) p->poly_vert[i] = GEO_Make_Geo_Coord(2);
1677 return(p);
1678 }
1679
1680 //////////////////////////////////////////////////////////////
1681 //////////////////////////////////////////////////////////////
1682
Make_Sarea(int n_poly)1683 t_sarea *Make_Sarea(int n_poly)
1684 {
1685 t_sarea *s;
1686 s = (t_sarea *)mCalloc(1,sizeof(t_sarea ));
1687 s->a_poly = (t_poly **)mCalloc(n_poly,sizeof(t_poly *));
1688 return(s);
1689 }
1690
1691 //////////////////////////////////////////////////////////////
1692 //////////////////////////////////////////////////////////////
1693
Make_Linked_List()1694 t_ll *Make_Linked_List()
1695 {
1696 t_ll *list;
1697
1698 // Create list if non-existing and add
1699 list = (t_ll *)mCalloc(1,sizeof(t_ll));
1700
1701 return list;
1702 }
1703
1704 //////////////////////////////////////////////////////////////
1705 //////////////////////////////////////////////////////////////
1706
1707 /* Matrices used in transfer bootstrap computation (tbe.c) */
Alloc_TBE_Matrices(int n_otu,short unsigned *** i_matrix,short unsigned *** c_matrix,short unsigned *** hamming,short unsigned ** min_dist,short unsigned ** min_dist_edge,int ** cluster_sizes)1708 void Alloc_TBE_Matrices(int n_otu, short unsigned*** i_matrix, short unsigned*** c_matrix,short unsigned*** hamming, short unsigned** min_dist, short unsigned** min_dist_edge, int** cluster_sizes){
1709 int i;
1710 int nb_edges = 2*n_otu-3;
1711 (*min_dist) = (short unsigned*) malloc(nb_edges*sizeof(short unsigned)); /* array of min Hamming distances */
1712 (*min_dist_edge) = (short unsigned*) malloc(nb_edges*sizeof(short unsigned)); /* array of edge ids corresponding to min Hamming distances */
1713 (*cluster_sizes) = (int*) malloc(nb_edges*sizeof(int)); /* array of sizes of clusters associated to each branch (in the post order traversal) */
1714 (*c_matrix) = (short unsigned**) malloc(nb_edges*sizeof(short unsigned*)); /* matrix of cardinals of complements */
1715 (*i_matrix) = (short unsigned**) malloc(nb_edges*sizeof(short unsigned*)); /* matrix of cardinals of intersections */
1716 (*hamming) = (short unsigned**) malloc(nb_edges*sizeof(short unsigned*)); /* matrix of Hamming distances */
1717 for (i=0; i<nb_edges; i++){
1718 (*c_matrix)[i] = (short unsigned*) malloc(nb_edges*sizeof(short unsigned));
1719 (*i_matrix)[i] = (short unsigned*) malloc(nb_edges*sizeof(short unsigned));
1720 (*hamming)[i] = (short unsigned*) malloc(nb_edges*sizeof(short unsigned));
1721 (*min_dist)[i] = n_otu; /* initialization to the nb of taxa */
1722 (*cluster_sizes)[i] = 0;
1723 }
1724 }
1725
1726 //////////////////////////////////////////////////////////////
1727 //////////////////////////////////////////////////////////////
1728
Make_MutMap(t_tree * tree)1729 void Make_MutMap(t_tree *tree)
1730 {
1731 // (# of edges) X (# of sites) X (# mutation types: A<->C A<->G A<->T C<->G C<->T G<->T)
1732 do
1733 {
1734 tree->mutmap = (int *)mCalloc((2*tree->n_otu-3)*(tree->n_pattern)*6,sizeof(int));
1735 tree = tree->next_mixt;
1736 }
1737 while(tree);
1738 }
1739
1740 //////////////////////////////////////////////////////////////
1741 //////////////////////////////////////////////////////////////
1742
Make_Label()1743 t_label *Make_Label()
1744 {
1745 t_label *lab;
1746
1747 lab = (t_label *)mCalloc(1,sizeof(t_label));
1748 lab->key = (char *)mCalloc(T_MAX_KEY,sizeof(char));
1749 lab->val = (char *)mCalloc(T_MAX_VAL,sizeof(char));
1750 lab->next = NULL;
1751
1752 return(lab);
1753 }
1754
1755 //////////////////////////////////////////////////////////////
1756 //////////////////////////////////////////////////////////////
1757
Make_Contrasts(t_tree * tree)1758 void Make_Contrasts(t_tree *tree)
1759 {
1760 tree->ctrst = (t_ctrst *)mCalloc(1,sizeof(t_ctrst));
1761 tree->ctrst->x = (phydbl *)mCalloc(2*tree->n_otu,sizeof(phydbl));
1762 tree->ctrst->tprime = (phydbl *)mCalloc(2*tree->n_otu,sizeof(phydbl));
1763 }
1764
1765 //////////////////////////////////////////////////////////////
1766 //////////////////////////////////////////////////////////////
1767
1768 //////////////////////////////////////////////////////////////
1769 //////////////////////////////////////////////////////////////
1770
1771 //////////////////////////////////////////////////////////////
1772 //////////////////////////////////////////////////////////////
1773
1774 //////////////////////////////////////////////////////////////
1775 //////////////////////////////////////////////////////////////
1776
1777 //////////////////////////////////////////////////////////////
1778 //////////////////////////////////////////////////////////////
1779
1780 //////////////////////////////////////////////////////////////
1781 //////////////////////////////////////////////////////////////
1782
1783 //////////////////////////////////////////////////////////////
1784 //////////////////////////////////////////////////////////////
1785
1786 //////////////////////////////////////////////////////////////
1787 //////////////////////////////////////////////////////////////
1788
1789 //////////////////////////////////////////////////////////////
1790 //////////////////////////////////////////////////////////////
1791
1792 //////////////////////////////////////////////////////////////
1793 //////////////////////////////////////////////////////////////
1794
1795 //////////////////////////////////////////////////////////////
1796 //////////////////////////////////////////////////////////////
1797
1798 //////////////////////////////////////////////////////////////
1799 //////////////////////////////////////////////////////////////
1800
1801 //////////////////////////////////////////////////////////////
1802 //////////////////////////////////////////////////////////////
1803
1804
1805
1806