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 | |
---|
17 | void Make_Tree_4_Lk(t_tree *tree, calign *cdata, int n_site) |
---|
18 | { |
---|
19 | int i; |
---|
20 | |
---|
21 | tree->c_lnL_sorted = (phydbl *)mCalloc(tree->n_pattern, sizeof(phydbl)); |
---|
22 | tree->cur_site_lk = (phydbl *)mCalloc(cdata->crunch_len,sizeof(phydbl)); |
---|
23 | tree->old_site_lk = (phydbl *)mCalloc(cdata->crunch_len,sizeof(phydbl)); |
---|
24 | tree->site_lk_cat = (phydbl *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(phydbl)); |
---|
25 | tree->log_site_lk_cat = (phydbl **)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(phydbl *)); |
---|
26 | For(i,MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes)) |
---|
27 | tree->log_site_lk_cat[i] = (phydbl *)mCalloc(cdata->crunch_len,sizeof(phydbl)); |
---|
28 | |
---|
29 | tree->log_lks_aLRT = (phydbl **)mCalloc(3,sizeof(phydbl *)); |
---|
30 | For(i,3) tree->log_lks_aLRT[i] = (phydbl *)mCalloc(tree->data->init_len,sizeof(phydbl)); |
---|
31 | |
---|
32 | For(i,2*tree->n_otu-1) Make_Edge_NNI(tree->a_edges[i]); |
---|
33 | |
---|
34 | if(tree->is_mixt_tree == NO) |
---|
35 | { |
---|
36 | For(i,2*tree->n_otu-1) Make_Edge_Lk(tree->a_edges[i],tree); |
---|
37 | For(i,2*tree->n_otu-2) Make_Node_Lk(tree->a_nodes[i]); |
---|
38 | |
---|
39 | if(tree->mod->s_opt->greedy) |
---|
40 | Init_P_Lk_Tips_Double(tree); |
---|
41 | else |
---|
42 | Init_P_Lk_Tips_Int(tree); |
---|
43 | |
---|
44 | Init_P_Lk_Loc(tree); |
---|
45 | } |
---|
46 | } |
---|
47 | |
---|
48 | ////////////////////////////////////////////////////////////// |
---|
49 | ////////////////////////////////////////////////////////////// |
---|
50 | |
---|
51 | |
---|
52 | void Make_Tree_4_Pars(t_tree *tree, calign *cdata, int n_site) |
---|
53 | { |
---|
54 | int i; |
---|
55 | tree->site_pars = (int *)mCalloc(tree->n_pattern,sizeof(int)); |
---|
56 | tree->step_mat = (int *)mCalloc(tree->mod->ns * tree->mod->ns,sizeof(int)); |
---|
57 | For(i,2*tree->n_otu-1) Make_Edge_Pars(tree->a_edges[i],tree); |
---|
58 | Init_Ui_Tips(tree); |
---|
59 | Init_P_Pars_Tips(tree); /* Must be called after Init_Ui_Tips is called */ |
---|
60 | Get_Step_Mat(tree); |
---|
61 | } |
---|
62 | |
---|
63 | ////////////////////////////////////////////////////////////// |
---|
64 | ////////////////////////////////////////////////////////////// |
---|
65 | |
---|
66 | void Make_All_Edges_Lk(t_node *a, t_node *d, t_tree *tree) |
---|
67 | { |
---|
68 | int i; |
---|
69 | |
---|
70 | For(i,3) if((a->v[i]) && (a->v[i] == d)) Make_Edge_Lk(a->b[i],tree); |
---|
71 | if(d->tax) return; |
---|
72 | else |
---|
73 | { |
---|
74 | For(i,3) |
---|
75 | { |
---|
76 | if(d->v[i] != a) |
---|
77 | Make_All_Edges_Lk(d,d->v[i],tree); |
---|
78 | } |
---|
79 | } |
---|
80 | } |
---|
81 | |
---|
82 | ////////////////////////////////////////////////////////////// |
---|
83 | ////////////////////////////////////////////////////////////// |
---|
84 | |
---|
85 | void Make_New_Edge_Label(t_edge *b) |
---|
86 | { |
---|
87 | int i; |
---|
88 | |
---|
89 | b->labels = (char **)realloc(b->labels,(b->n_labels+BLOCK_LABELS)*sizeof(char *)); |
---|
90 | |
---|
91 | if(!b->labels) |
---|
92 | { |
---|
93 | PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__); |
---|
94 | Warn_And_Exit(""); |
---|
95 | } |
---|
96 | else |
---|
97 | { |
---|
98 | for(i=b->n_labels;i<b->n_labels+BLOCK_LABELS;i++) b->labels[i] = (char *)mCalloc(T_MAX_LABEL,sizeof(char)); |
---|
99 | } |
---|
100 | } |
---|
101 | |
---|
102 | ////////////////////////////////////////////////////////////// |
---|
103 | ////////////////////////////////////////////////////////////// |
---|
104 | |
---|
105 | t_edge *Make_Edge_Light(t_node *a, t_node *d, int num) |
---|
106 | { |
---|
107 | t_edge *b; |
---|
108 | |
---|
109 | b = (t_edge *)mCalloc(1,sizeof(t_edge)); |
---|
110 | |
---|
111 | b->l = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
112 | Init_Scalar_Dbl(b->l); |
---|
113 | |
---|
114 | b->l_old = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
115 | Init_Scalar_Dbl(b->l_old); |
---|
116 | |
---|
117 | b->l_var = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
118 | Init_Scalar_Dbl(b->l_var); |
---|
119 | |
---|
120 | b->l_var_old = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
121 | Init_Scalar_Dbl(b->l_var_old); |
---|
122 | |
---|
123 | Init_Edge_Light(b,num); |
---|
124 | |
---|
125 | if(a && b) |
---|
126 | { |
---|
127 | b->left = a; b->rght = d; |
---|
128 | if(a->tax) {b->rght = a; b->left = d;} /* root */ |
---|
129 | /* a tip is necessary on the right side of the t_edge */ |
---|
130 | |
---|
131 | (b->left == a)? |
---|
132 | (Make_Edge_Dirs(b,a,d,NULL)): |
---|
133 | (Make_Edge_Dirs(b,d,a,NULL)); |
---|
134 | |
---|
135 | b->l->v = a->l[b->l_r]; |
---|
136 | if(a->tax) b->l->v = a->l[b->r_l]; |
---|
137 | b->l_old->v = b->l->v; |
---|
138 | } |
---|
139 | else |
---|
140 | { |
---|
141 | b->left = NULL; |
---|
142 | b->rght = NULL; |
---|
143 | } |
---|
144 | |
---|
145 | return b; |
---|
146 | |
---|
147 | } |
---|
148 | |
---|
149 | ////////////////////////////////////////////////////////////// |
---|
150 | ////////////////////////////////////////////////////////////// |
---|
151 | |
---|
152 | void Make_Edge_Pars(t_edge *b, t_tree *tree) |
---|
153 | { |
---|
154 | Make_Edge_Pars_Left(b,tree); |
---|
155 | Make_Edge_Pars_Rght(b,tree); |
---|
156 | } |
---|
157 | |
---|
158 | ////////////////////////////////////////////////////////////// |
---|
159 | ////////////////////////////////////////////////////////////// |
---|
160 | |
---|
161 | void Make_Edge_Pars_Left(t_edge *b, t_tree *tree) |
---|
162 | { |
---|
163 | b->pars_l = (int *)mCalloc(tree->data->crunch_len,sizeof(int)); |
---|
164 | b->ui_l = (unsigned int *)mCalloc(tree->data->crunch_len,sizeof(unsigned int)); |
---|
165 | b->p_pars_l = (int *)mCalloc(tree->data->crunch_len*tree->mod->ns,sizeof(int )); |
---|
166 | } |
---|
167 | |
---|
168 | ////////////////////////////////////////////////////////////// |
---|
169 | ////////////////////////////////////////////////////////////// |
---|
170 | |
---|
171 | void Make_Edge_Pars_Rght(t_edge *b, t_tree *tree) |
---|
172 | { |
---|
173 | b->pars_r = (int *)mCalloc(tree->data->crunch_len,sizeof(int)); |
---|
174 | b->ui_r = (unsigned int *)mCalloc(tree->data->crunch_len,sizeof(unsigned int)); |
---|
175 | b->p_pars_r = (int *)mCalloc(tree->data->crunch_len*tree->mod->ns,sizeof(int )); |
---|
176 | } |
---|
177 | |
---|
178 | ////////////////////////////////////////////////////////////// |
---|
179 | ////////////////////////////////////////////////////////////// |
---|
180 | |
---|
181 | void Make_Edge_Lk(t_edge *b, t_tree *tree) |
---|
182 | { |
---|
183 | if(tree->is_mixt_tree) |
---|
184 | { |
---|
185 | PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__); |
---|
186 | Warn_And_Exit(""); |
---|
187 | } |
---|
188 | |
---|
189 | b->l_old->v = b->l->v; |
---|
190 | |
---|
191 | b->Pij_rr = (phydbl *)mCalloc(tree->mod->ras->n_catg*tree->mod->ns*tree->mod->ns,sizeof(phydbl)); |
---|
192 | |
---|
193 | Make_Edge_Lk_Left(b,tree); |
---|
194 | Make_Edge_Lk_Rght(b,tree); |
---|
195 | } |
---|
196 | |
---|
197 | ////////////////////////////////////////////////////////////// |
---|
198 | ////////////////////////////////////////////////////////////// |
---|
199 | |
---|
200 | void Make_Edge_Lk_Left(t_edge *b, t_tree *tree) |
---|
201 | { |
---|
202 | int ns = tree->mod->ns; |
---|
203 | |
---|
204 | b->div_post_pred_left = (short int *)mCalloc(ns,sizeof(short int)); |
---|
205 | |
---|
206 | b->sum_scale_left_cat = (int *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int)); |
---|
207 | |
---|
208 | if(b->left && !b->left->tax) |
---|
209 | b->sum_scale_left = (int *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int)); |
---|
210 | else |
---|
211 | b->sum_scale_left = NULL; |
---|
212 | |
---|
213 | |
---|
214 | if(b->left) |
---|
215 | { |
---|
216 | if((!b->left->tax) || (tree->mod->s_opt->greedy)) |
---|
217 | { |
---|
218 | 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)); |
---|
219 | b->p_lk_tip_l = NULL; |
---|
220 | } |
---|
221 | else if(b->left->tax) |
---|
222 | { |
---|
223 | b->p_lk_left = NULL; |
---|
224 | b->p_lk_tip_l = (short int *)mCalloc(tree->data->crunch_len*tree->mod->ns,sizeof(short int )); |
---|
225 | } |
---|
226 | } |
---|
227 | else |
---|
228 | { |
---|
229 | b->p_lk_left = NULL; |
---|
230 | b->p_lk_tip_l = NULL; |
---|
231 | } |
---|
232 | |
---|
233 | if(b->num >= 2*tree->n_otu-3) |
---|
234 | { |
---|
235 | b->sum_scale_left = (int *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int)); |
---|
236 | 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)); |
---|
237 | } |
---|
238 | |
---|
239 | |
---|
240 | b->patt_id_left = (int *)mCalloc(tree->data->crunch_len,sizeof(int)); |
---|
241 | b->p_lk_loc_left = (int *)mCalloc(tree->data->crunch_len,sizeof(int)); |
---|
242 | } |
---|
243 | |
---|
244 | ////////////////////////////////////////////////////////////// |
---|
245 | ////////////////////////////////////////////////////////////// |
---|
246 | |
---|
247 | void Make_Edge_Lk_Rght(t_edge *b, t_tree *tree) |
---|
248 | { |
---|
249 | int ns = tree->mod->ns; |
---|
250 | |
---|
251 | b->div_post_pred_rght = (short int *)mCalloc(ns,sizeof(short int)); |
---|
252 | |
---|
253 | b->sum_scale_rght_cat = (int *)mCalloc(MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int)); |
---|
254 | |
---|
255 | if(b->rght && !b->rght->tax) |
---|
256 | b->sum_scale_rght = (int *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int)); |
---|
257 | else |
---|
258 | b->sum_scale_rght = NULL; |
---|
259 | |
---|
260 | |
---|
261 | if(b->rght) |
---|
262 | { |
---|
263 | if((!b->rght->tax) || (tree->mod->s_opt->greedy)) |
---|
264 | { |
---|
265 | 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)); |
---|
266 | b->p_lk_tip_r = NULL; |
---|
267 | } |
---|
268 | else if(b->rght->tax) |
---|
269 | { |
---|
270 | b->p_lk_rght = NULL; |
---|
271 | b->p_lk_tip_r = (short int *)mCalloc(tree->data->crunch_len*tree->mod->ns,sizeof(short int )); |
---|
272 | } |
---|
273 | } |
---|
274 | else |
---|
275 | { |
---|
276 | b->p_lk_rght = NULL; |
---|
277 | b->p_lk_tip_r = NULL; |
---|
278 | } |
---|
279 | |
---|
280 | if(b->num >= 2*tree->n_otu-3) |
---|
281 | { |
---|
282 | b->sum_scale_rght = (int *)mCalloc(tree->data->crunch_len*MAX(tree->mod->ras->n_catg,tree->mod->n_mixt_classes),sizeof(int)); |
---|
283 | 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)); |
---|
284 | } |
---|
285 | |
---|
286 | b->patt_id_rght = (int *)mCalloc(tree->data->crunch_len,sizeof(int)); |
---|
287 | b->p_lk_loc_rght = (int *)mCalloc(tree->data->crunch_len,sizeof(int)); |
---|
288 | } |
---|
289 | |
---|
290 | ////////////////////////////////////////////////////////////// |
---|
291 | ////////////////////////////////////////////////////////////// |
---|
292 | |
---|
293 | void Make_Edge_NNI(t_edge *b) |
---|
294 | { |
---|
295 | b->nni = Make_NNI(); |
---|
296 | b->nni->b = b; |
---|
297 | b->nni->left = b->left; |
---|
298 | b->nni->rght = b->rght; |
---|
299 | } |
---|
300 | |
---|
301 | ////////////////////////////////////////////////////////////// |
---|
302 | ////////////////////////////////////////////////////////////// |
---|
303 | |
---|
304 | |
---|
305 | nni *Make_NNI() |
---|
306 | { |
---|
307 | nni *a_nni; |
---|
308 | a_nni = (nni *)mCalloc(1,sizeof(nni )); |
---|
309 | Init_NNI(a_nni); |
---|
310 | return a_nni; |
---|
311 | } |
---|
312 | |
---|
313 | ////////////////////////////////////////////////////////////// |
---|
314 | ////////////////////////////////////////////////////////////// |
---|
315 | |
---|
316 | t_node *Make_Node_Light(int num) |
---|
317 | { |
---|
318 | t_node *n; |
---|
319 | |
---|
320 | n = (t_node *)mCalloc(1,sizeof(t_node)); |
---|
321 | n->v = (t_node **)mCalloc(3,sizeof(t_node *)); |
---|
322 | n->l = (phydbl *)mCalloc(3,sizeof(phydbl)); |
---|
323 | n->b = (t_edge **)mCalloc(3,sizeof(t_edge *)); |
---|
324 | n->score = (phydbl *)mCalloc(3,sizeof(phydbl)); |
---|
325 | n->s_ingrp = (int *)mCalloc(3,sizeof(int)); |
---|
326 | n->s_outgrp = (int *)mCalloc(3,sizeof(int)); |
---|
327 | |
---|
328 | Init_Node_Light(n,num); |
---|
329 | |
---|
330 | return n; |
---|
331 | } |
---|
332 | |
---|
333 | ////////////////////////////////////////////////////////////// |
---|
334 | ////////////////////////////////////////////////////////////// |
---|
335 | |
---|
336 | |
---|
337 | |
---|
338 | void Make_Node_Lk(t_node *n) |
---|
339 | { |
---|
340 | /* n->n_ex_nodes = (int *)mCalloc(2,sizeof(int)); */ |
---|
341 | return; |
---|
342 | } |
---|
343 | |
---|
344 | nexcom **Make_Nexus_Com() |
---|
345 | { |
---|
346 | nexcom **com; |
---|
347 | int i; |
---|
348 | |
---|
349 | com = (nexcom **)mCalloc(N_MAX_NEX_COM,sizeof(nexcom *)); |
---|
350 | |
---|
351 | For(i,N_MAX_NEX_COM) |
---|
352 | { |
---|
353 | com[i] = (nexcom *)mCalloc(1,sizeof(nexcom)); |
---|
354 | com[i]->name = (char *)mCalloc(T_MAX_NEX_COM,sizeof(char)); |
---|
355 | com[i]->parm = (nexparm **)mCalloc(N_MAX_NEX_PARM,sizeof(nexparm *)); |
---|
356 | } |
---|
357 | |
---|
358 | return com; |
---|
359 | } |
---|
360 | |
---|
361 | ////////////////////////////////////////////////////////////// |
---|
362 | ////////////////////////////////////////////////////////////// |
---|
363 | |
---|
364 | |
---|
365 | nexparm *Make_Nexus_Parm() |
---|
366 | { |
---|
367 | nexparm *parm; |
---|
368 | |
---|
369 | parm = (nexparm *)mCalloc(1,sizeof(nexparm)); |
---|
370 | parm->name = (char *)mCalloc(T_MAX_TOKEN,sizeof(char )); |
---|
371 | parm->value = (char *)mCalloc(T_MAX_TOKEN,sizeof(char )); |
---|
372 | |
---|
373 | return parm; |
---|
374 | } |
---|
375 | |
---|
376 | ////////////////////////////////////////////////////////////// |
---|
377 | ////////////////////////////////////////////////////////////// |
---|
378 | |
---|
379 | matrix *Make_Mat(int n_otu) |
---|
380 | { |
---|
381 | matrix *mat; |
---|
382 | int i; |
---|
383 | |
---|
384 | mat = (matrix *)mCalloc(1,sizeof(matrix)); |
---|
385 | |
---|
386 | mat->n_otu = n_otu; |
---|
387 | |
---|
388 | mat->P = (phydbl **)mCalloc(n_otu,sizeof(phydbl *)); |
---|
389 | mat->Q = (phydbl **)mCalloc(n_otu,sizeof(phydbl *)); |
---|
390 | mat->dist = (phydbl **)mCalloc(n_otu,sizeof(phydbl *)); |
---|
391 | mat->on_off = (int *)mCalloc(n_otu,sizeof(int)); |
---|
392 | mat->name = (char **)mCalloc(n_otu,sizeof(char *)); |
---|
393 | mat->tip_node = (t_node **)mCalloc(n_otu,sizeof(t_node *)); |
---|
394 | |
---|
395 | |
---|
396 | For(i,n_otu) |
---|
397 | { |
---|
398 | mat->P[i] = (phydbl *)mCalloc(n_otu,sizeof(phydbl)); |
---|
399 | mat->Q[i] = (phydbl *)mCalloc(n_otu,sizeof(phydbl)); |
---|
400 | mat->dist[i] = (phydbl *)mCalloc(n_otu,sizeof(phydbl)); |
---|
401 | mat->name[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char)); |
---|
402 | } |
---|
403 | |
---|
404 | return mat; |
---|
405 | } |
---|
406 | |
---|
407 | ////////////////////////////////////////////////////////////// |
---|
408 | ////////////////////////////////////////////////////////////// |
---|
409 | |
---|
410 | t_tree *Make_Tree_From_Scratch(int n_otu, calign *data) |
---|
411 | { |
---|
412 | t_tree *tree; |
---|
413 | |
---|
414 | tree = Make_Tree(n_otu); |
---|
415 | Init_Tree(tree,n_otu); |
---|
416 | Make_All_Tree_Nodes(tree); |
---|
417 | Make_All_Tree_Edges(tree); |
---|
418 | Make_Tree_Path(tree); |
---|
419 | if(data) |
---|
420 | { |
---|
421 | Copy_Tax_Names_To_Tip_Labels(tree,data); |
---|
422 | tree->data = data; |
---|
423 | } |
---|
424 | return tree; |
---|
425 | } |
---|
426 | |
---|
427 | ////////////////////////////////////////////////////////////// |
---|
428 | ////////////////////////////////////////////////////////////// |
---|
429 | |
---|
430 | |
---|
431 | t_tree *Make_Tree(int n_otu) |
---|
432 | { |
---|
433 | t_tree *tree; |
---|
434 | tree = (t_tree *)mCalloc(1,sizeof(t_tree )); |
---|
435 | tree->t_dir = (short int *)mCalloc((2*n_otu-2)*(2*n_otu-2),sizeof(int)); |
---|
436 | return tree; |
---|
437 | } |
---|
438 | |
---|
439 | ////////////////////////////////////////////////////////////// |
---|
440 | ////////////////////////////////////////////////////////////// |
---|
441 | |
---|
442 | |
---|
443 | void Make_Tree_Path(t_tree *tree) |
---|
444 | { |
---|
445 | tree->curr_path = (t_node **)mCalloc(tree->n_otu,sizeof(t_node *)); |
---|
446 | } |
---|
447 | |
---|
448 | ////////////////////////////////////////////////////////////// |
---|
449 | ////////////////////////////////////////////////////////////// |
---|
450 | |
---|
451 | |
---|
452 | void Make_All_Tree_Nodes(t_tree *tree) |
---|
453 | { |
---|
454 | int i; |
---|
455 | |
---|
456 | tree->a_nodes = (t_node **)mCalloc(2*tree->n_otu-1,sizeof(t_node *)); |
---|
457 | |
---|
458 | For(i,2*tree->n_otu-1) |
---|
459 | { |
---|
460 | tree->a_nodes[i] = (t_node *)Make_Node_Light(i); |
---|
461 | if(i < tree->n_otu) tree->a_nodes[i]->tax = 1; |
---|
462 | else tree->a_nodes[i]->tax = 0; |
---|
463 | } |
---|
464 | } |
---|
465 | |
---|
466 | ////////////////////////////////////////////////////////////// |
---|
467 | ////////////////////////////////////////////////////////////// |
---|
468 | |
---|
469 | void Make_All_Tree_Edges(t_tree *tree) |
---|
470 | { |
---|
471 | int i; |
---|
472 | tree->a_edges = (t_edge **)mCalloc(2*tree->n_otu-1,sizeof(t_edge *)); |
---|
473 | For(i,2*tree->n_otu-1) tree->a_edges[i] = (t_edge *)Make_Edge_Light(NULL,NULL,i); |
---|
474 | } |
---|
475 | |
---|
476 | ////////////////////////////////////////////////////////////// |
---|
477 | ////////////////////////////////////////////////////////////// |
---|
478 | |
---|
479 | calign *Make_Cseq(int n_otu, int crunch_len, int state_len, int init_len, char **sp_names) |
---|
480 | { |
---|
481 | calign *cdata; |
---|
482 | int j; |
---|
483 | |
---|
484 | cdata = (calign *)mCalloc(1,sizeof(calign)); |
---|
485 | cdata->n_otu = n_otu; |
---|
486 | cdata->c_seq = (align **)mCalloc(n_otu,sizeof(align *)); |
---|
487 | cdata->b_frq = (phydbl *)mCalloc(T_MAX_ALPHABET,sizeof(phydbl)); |
---|
488 | cdata->wght = (int *)mCalloc(crunch_len,sizeof(int)); |
---|
489 | cdata->ambigu = (short int *)mCalloc(crunch_len,sizeof(short int)); |
---|
490 | cdata->invar = (short int *)mCalloc(crunch_len,sizeof(short int)); |
---|
491 | cdata->sitepatt = (int *)mCalloc(init_len,sizeof(int )); |
---|
492 | cdata->format = 0; |
---|
493 | |
---|
494 | cdata->crunch_len = crunch_len; |
---|
495 | cdata->init_len = init_len; |
---|
496 | cdata->obs_pinvar = .0; |
---|
497 | |
---|
498 | For(j,n_otu) |
---|
499 | { |
---|
500 | cdata->c_seq[j] = (align *)mCalloc(1,sizeof(align)); |
---|
501 | cdata->c_seq[j]->name = (char *)mCalloc((int)(strlen(sp_names[j])+1),sizeof(char)); |
---|
502 | strcpy(cdata->c_seq[j]->name,sp_names[j]); |
---|
503 | cdata->c_seq[j]->state = (char *)mCalloc(crunch_len*state_len,sizeof(char)); |
---|
504 | cdata->c_seq[j]->d_state = (int *)mCalloc(crunch_len*state_len,sizeof(int)); |
---|
505 | cdata->c_seq[j]->is_ambigu = (short int *)mCalloc(crunch_len,sizeof(short int)); |
---|
506 | } |
---|
507 | |
---|
508 | return cdata; |
---|
509 | } |
---|
510 | |
---|
511 | ////////////////////////////////////////////////////////////// |
---|
512 | ////////////////////////////////////////////////////////////// |
---|
513 | |
---|
514 | |
---|
515 | t_treelist *Make_Treelist(int list_size) |
---|
516 | { |
---|
517 | t_treelist *tlist; |
---|
518 | |
---|
519 | tlist = (t_treelist *)mCalloc(1,sizeof(t_treelist)); |
---|
520 | tlist->list_size = list_size; |
---|
521 | tlist->tree = (t_tree **)mCalloc(list_size,sizeof(t_tree *)); |
---|
522 | |
---|
523 | return tlist; |
---|
524 | } |
---|
525 | |
---|
526 | ////////////////////////////////////////////////////////////// |
---|
527 | ////////////////////////////////////////////////////////////// |
---|
528 | |
---|
529 | t_opt *Make_Optimiz() |
---|
530 | { |
---|
531 | t_opt *s_opt; |
---|
532 | s_opt = (t_opt *)mCalloc(1,sizeof(t_opt)); |
---|
533 | return s_opt; |
---|
534 | } |
---|
535 | |
---|
536 | ////////////////////////////////////////////////////////////// |
---|
537 | ////////////////////////////////////////////////////////////// |
---|
538 | |
---|
539 | void Make_Custom_Model(t_mod *mod) |
---|
540 | { |
---|
541 | if(!mod->r_mat) |
---|
542 | { |
---|
543 | PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__); |
---|
544 | Warn_And_Exit(""); |
---|
545 | } |
---|
546 | |
---|
547 | if(!mod->r_mat->rr->v) |
---|
548 | mod->r_mat->rr->v = (phydbl *)mCalloc(mod->ns*(mod->ns-1)/2,sizeof(phydbl)); |
---|
549 | |
---|
550 | if(!mod->r_mat->rr_val->v) |
---|
551 | mod->r_mat->rr_val->v = (phydbl *)mCalloc(mod->ns*(mod->ns-1)/2,sizeof(phydbl)); |
---|
552 | |
---|
553 | if(!mod->r_mat->rr_num->v) |
---|
554 | mod->r_mat->rr_num->v = (int *)mCalloc(mod->ns*(mod->ns-1)/2,sizeof(int *)); |
---|
555 | |
---|
556 | if(!mod->r_mat->n_rr_per_cat->v) |
---|
557 | mod->r_mat->n_rr_per_cat->v = (int *)mCalloc(mod->ns*(mod->ns-1)/2,sizeof(int)); |
---|
558 | } |
---|
559 | |
---|
560 | ////////////////////////////////////////////////////////////// |
---|
561 | ////////////////////////////////////////////////////////////// |
---|
562 | |
---|
563 | t_string *Make_String(int len) |
---|
564 | { |
---|
565 | t_string *ts; |
---|
566 | |
---|
567 | ts = (t_string *)mCalloc(1,sizeof(t_string)); |
---|
568 | ts->s = (char *)mCalloc(len,sizeof(char)); |
---|
569 | |
---|
570 | return(ts); |
---|
571 | } |
---|
572 | |
---|
573 | ////////////////////////////////////////////////////////////// |
---|
574 | ////////////////////////////////////////////////////////////// |
---|
575 | |
---|
576 | t_mod *Make_Model_Basic() |
---|
577 | { |
---|
578 | t_mod *mod; |
---|
579 | |
---|
580 | mod = (t_mod *)mCalloc(1,sizeof(t_mod)); |
---|
581 | |
---|
582 | mod->modelname = Make_String(T_MAX_NAME); |
---|
583 | Init_String(mod->modelname); |
---|
584 | |
---|
585 | mod->custom_mod_string = Make_String(T_MAX_NAME); |
---|
586 | Init_String(mod->custom_mod_string); |
---|
587 | |
---|
588 | mod->ras = Make_RAS_Basic(); |
---|
589 | |
---|
590 | mod->kappa = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
591 | Init_Scalar_Dbl(mod->kappa); |
---|
592 | |
---|
593 | mod->lambda = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
594 | Init_Scalar_Dbl(mod->lambda); |
---|
595 | |
---|
596 | mod->br_len_multiplier = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
597 | Init_Scalar_Dbl(mod->br_len_multiplier); |
---|
598 | |
---|
599 | mod->mr = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
600 | Init_Scalar_Dbl(mod->mr); |
---|
601 | |
---|
602 | mod->user_b_freq = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
603 | Init_Vect_Dbl(0,mod->user_b_freq); |
---|
604 | mod->user_b_freq->v = (phydbl *)mCalloc(T_MAX_OPTION,sizeof(phydbl)); |
---|
605 | |
---|
606 | mod->e_frq_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
607 | Init_Scalar_Dbl(mod->e_frq_weight); |
---|
608 | |
---|
609 | mod->r_mat_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
610 | Init_Scalar_Dbl(mod->r_mat_weight); |
---|
611 | |
---|
612 | mod->aa_rate_mat_file = Make_String(T_MAX_FILE); |
---|
613 | Init_String(mod->aa_rate_mat_file); |
---|
614 | return mod; |
---|
615 | } |
---|
616 | |
---|
617 | ////////////////////////////////////////////////////////////// |
---|
618 | ////////////////////////////////////////////////////////////// |
---|
619 | |
---|
620 | /*! Call only when the values of mod->ns & ras->n_catg is set to its final value */ |
---|
621 | |
---|
622 | void Make_Model_Complete(t_mod *mod) |
---|
623 | { |
---|
624 | if(mod->use_m4mod == YES) |
---|
625 | { |
---|
626 | M4_Make_Complete(mod->m4mod->n_h,mod->m4mod->n_o,mod->m4mod); |
---|
627 | mod->ns = mod->m4mod->n_o * mod->m4mod->n_h; |
---|
628 | } |
---|
629 | |
---|
630 | mod->Pij_rr = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
631 | Init_Vect_Dbl(0,mod->Pij_rr); |
---|
632 | mod->Pij_rr->v = (phydbl *)mCalloc(mod->ras->n_catg*mod->ns*mod->ns,sizeof(phydbl)); |
---|
633 | |
---|
634 | mod->eigen = (eigen *)Make_Eigen_Struct(mod->ns); |
---|
635 | |
---|
636 | // If r_mat (e_frq) are not NULL, then they have been created elsewhere and affected. |
---|
637 | if(!mod->r_mat) |
---|
638 | { |
---|
639 | mod->r_mat = (t_rmat *)Make_Rmat(mod->ns); |
---|
640 | Init_Rmat(mod->r_mat); |
---|
641 | } |
---|
642 | if(!mod->e_frq) |
---|
643 | { |
---|
644 | mod->e_frq = (t_efrq *)Make_Efrq(mod->ns); |
---|
645 | Init_Efrq(mod->e_frq); |
---|
646 | } |
---|
647 | |
---|
648 | Make_RAS_Complete(mod->ras); |
---|
649 | |
---|
650 | mod->user_b_freq->len = mod->ns; |
---|
651 | |
---|
652 | if(mod->whichmodel < 0) |
---|
653 | { |
---|
654 | PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__); |
---|
655 | Exit("\n"); |
---|
656 | } |
---|
657 | |
---|
658 | if(mod->whichmodel == CUSTOM) |
---|
659 | { |
---|
660 | Make_Custom_Model(mod); |
---|
661 | Translate_Custom_Mod_String(mod); |
---|
662 | } |
---|
663 | if(mod->whichmodel == GTR) |
---|
664 | { |
---|
665 | Make_Custom_Model(mod); |
---|
666 | } |
---|
667 | } |
---|
668 | |
---|
669 | ////////////////////////////////////////////////////////////// |
---|
670 | ////////////////////////////////////////////////////////////// |
---|
671 | |
---|
672 | t_ras *Make_RAS_Basic() |
---|
673 | { |
---|
674 | t_ras *ras; |
---|
675 | |
---|
676 | ras = (t_ras *)mCalloc(1,sizeof(t_ras)); |
---|
677 | |
---|
678 | ras->gamma_r_proba = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
679 | Init_Vect_Dbl(0,ras->gamma_r_proba); |
---|
680 | ras->gamma_r_proba->v = NULL; |
---|
681 | |
---|
682 | ras->gamma_r_proba_unscaled = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
683 | Init_Vect_Dbl(0,ras->gamma_r_proba_unscaled); |
---|
684 | ras->gamma_r_proba_unscaled->v = NULL; |
---|
685 | |
---|
686 | ras->gamma_rr = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
687 | Init_Vect_Dbl(0,ras->gamma_rr); |
---|
688 | ras->gamma_rr->v = NULL; |
---|
689 | |
---|
690 | ras->gamma_rr_unscaled = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
691 | Init_Vect_Dbl(0,ras->gamma_rr_unscaled); |
---|
692 | ras->gamma_rr_unscaled->v = NULL; |
---|
693 | |
---|
694 | ras->alpha = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
695 | Init_Scalar_Dbl(ras->alpha); |
---|
696 | |
---|
697 | ras->pinvar = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
698 | Init_Scalar_Dbl(ras->pinvar); |
---|
699 | |
---|
700 | ras->free_rate_mr = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
701 | Init_Scalar_Dbl(ras->free_rate_mr); |
---|
702 | |
---|
703 | return(ras); |
---|
704 | } |
---|
705 | |
---|
706 | ////////////////////////////////////////////////////////////// |
---|
707 | ////////////////////////////////////////////////////////////// |
---|
708 | |
---|
709 | /*! Call only when the value of ras->n_catg is set to its final value */ |
---|
710 | void Make_RAS_Complete(t_ras *ras) |
---|
711 | { |
---|
712 | if(!ras->gamma_r_proba->v) |
---|
713 | { |
---|
714 | ras->gamma_r_proba->v = (phydbl *)mCalloc(ras->n_catg,sizeof(phydbl)); |
---|
715 | ras->gamma_r_proba_unscaled->v = (phydbl *)mCalloc(ras->n_catg,sizeof(phydbl)); |
---|
716 | ras->gamma_rr->v = (phydbl *)mCalloc(ras->n_catg,sizeof(phydbl)); |
---|
717 | ras->gamma_rr_unscaled->v = (phydbl *)mCalloc(ras->n_catg,sizeof(phydbl)); |
---|
718 | ras->skip_rate_cat = (short int *)mCalloc(ras->n_catg,sizeof(short int)); |
---|
719 | } |
---|
720 | } |
---|
721 | |
---|
722 | ////////////////////////////////////////////////////////////// |
---|
723 | ////////////////////////////////////////////////////////////// |
---|
724 | |
---|
725 | t_efrq *Make_Efrq(int ns) |
---|
726 | { |
---|
727 | t_efrq *e_frq; |
---|
728 | |
---|
729 | e_frq = (t_efrq *)mCalloc(1,sizeof(t_efrq)); |
---|
730 | |
---|
731 | e_frq->pi = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
732 | e_frq->pi->v = (phydbl *)mCalloc(ns,sizeof(phydbl)); |
---|
733 | e_frq->pi->len = ns; |
---|
734 | |
---|
735 | e_frq->pi_unscaled = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
736 | e_frq->pi_unscaled->v = (phydbl *)mCalloc(ns,sizeof(phydbl)); |
---|
737 | e_frq->pi_unscaled->len = ns; |
---|
738 | |
---|
739 | return e_frq; |
---|
740 | } |
---|
741 | |
---|
742 | ////////////////////////////////////////////////////////////// |
---|
743 | ////////////////////////////////////////////////////////////// |
---|
744 | |
---|
745 | t_rmat *Make_Rmat(int ns) |
---|
746 | { |
---|
747 | t_rmat *r_mat; |
---|
748 | |
---|
749 | r_mat = (t_rmat *)mCalloc(1,sizeof(t_rmat)); |
---|
750 | |
---|
751 | r_mat->qmat = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
752 | Init_Vect_Dbl(0,r_mat->qmat); |
---|
753 | |
---|
754 | r_mat->qmat_buff = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
755 | Init_Vect_Dbl(0,r_mat->qmat_buff); |
---|
756 | |
---|
757 | r_mat->rr = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
758 | Init_Vect_Dbl(0,r_mat->rr); |
---|
759 | |
---|
760 | r_mat->rr_val = (vect_dbl *)mCalloc(1,sizeof(vect_dbl)); |
---|
761 | Init_Vect_Dbl(0,r_mat->rr_val); |
---|
762 | |
---|
763 | r_mat->rr_num = (vect_int *)mCalloc(1,sizeof(vect_int)); |
---|
764 | Init_Vect_Int(0,r_mat->rr_num); |
---|
765 | |
---|
766 | r_mat->n_rr_per_cat = (vect_int *)mCalloc(1,sizeof(vect_int)); |
---|
767 | Init_Vect_Int(0,r_mat->n_rr_per_cat); |
---|
768 | |
---|
769 | r_mat->qmat->v = (phydbl *)mCalloc(ns*ns,sizeof(phydbl)); |
---|
770 | r_mat->qmat_buff->v = (phydbl *)mCalloc(ns*ns,sizeof(phydbl)); |
---|
771 | |
---|
772 | |
---|
773 | return(r_mat); |
---|
774 | } |
---|
775 | |
---|
776 | ////////////////////////////////////////////////////////////// |
---|
777 | ////////////////////////////////////////////////////////////// |
---|
778 | |
---|
779 | option *Make_Input() |
---|
780 | { |
---|
781 | int i; |
---|
782 | option* io = (option *)mCalloc(1,sizeof(option)); |
---|
783 | |
---|
784 | io->in_align_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
785 | io->in_tree_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
786 | io->in_constraint_tree_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
787 | io->out_tree_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
788 | io->out_trees_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
789 | io->out_boot_tree_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
790 | io->out_boot_stats_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
791 | io->out_stats_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
792 | io->out_lk_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
793 | io->out_ps_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
794 | io->out_trace_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
795 | io->nt_or_cd = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
796 | io->run_id_string = (char *)mCalloc(T_MAX_OPTION,sizeof(char)); |
---|
797 | io->clade_list_file = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
798 | io->alphabet = (char **)mCalloc(T_MAX_ALPHABET,sizeof(char *)); |
---|
799 | For(i,T_MAX_ALPHABET) io->alphabet[i] = (char *)mCalloc(T_MAX_STATE,sizeof(char )); |
---|
800 | io->treelist = (t_treelist *)mCalloc(1,sizeof(t_treelist)); |
---|
801 | io->mcmc = (t_mcmc *)MCMC_Make_MCMC_Struct(); |
---|
802 | io->rates = (t_rate *)RATES_Make_Rate_Struct(-1); |
---|
803 | |
---|
804 | return io; |
---|
805 | } |
---|
806 | |
---|
807 | ////////////////////////////////////////////////////////////// |
---|
808 | ////////////////////////////////////////////////////////////// |
---|
809 | |
---|
810 | t_mcmc *MCMC_Make_MCMC_Struct() |
---|
811 | { |
---|
812 | t_mcmc *mcmc; |
---|
813 | |
---|
814 | mcmc = (t_mcmc *)mCalloc(1,sizeof(t_mcmc)); |
---|
815 | mcmc->out_filename = (char *)mCalloc(T_MAX_FILE,sizeof(char)); |
---|
816 | |
---|
817 | return(mcmc); |
---|
818 | } |
---|
819 | |
---|
820 | ////////////////////////////////////////////////////////////// |
---|
821 | ////////////////////////////////////////////////////////////// |
---|
822 | |
---|
823 | eigen *Make_Eigen_Struct(int ns) |
---|
824 | { |
---|
825 | eigen *eig; |
---|
826 | |
---|
827 | eig = (eigen *)mCalloc(1,sizeof(eigen)); |
---|
828 | eig->size = ns; |
---|
829 | eig->space = (phydbl *)mCalloc(2*ns,sizeof(phydbl)); |
---|
830 | eig->space_int = (int *)mCalloc(2*ns,sizeof(int)); |
---|
831 | eig->e_val = (phydbl *)mCalloc(ns,sizeof(phydbl)); |
---|
832 | eig->e_val_im = (phydbl *)mCalloc(ns,sizeof(phydbl)); |
---|
833 | eig->r_e_vect = (phydbl *)mCalloc(ns*ns,sizeof(phydbl)); |
---|
834 | eig->r_e_vect_im = (phydbl *)mCalloc(ns*ns,sizeof(phydbl)); |
---|
835 | eig->l_e_vect = (phydbl *)mCalloc(ns*ns,sizeof(phydbl)); |
---|
836 | eig->q = (phydbl *)mCalloc(ns*ns,sizeof(phydbl)); |
---|
837 | |
---|
838 | Init_Eigen_Struct(eig); |
---|
839 | |
---|
840 | return eig; |
---|
841 | } |
---|
842 | |
---|
843 | ////////////////////////////////////////////////////////////// |
---|
844 | ////////////////////////////////////////////////////////////// |
---|
845 | |
---|
846 | triplet *Make_Triplet_Struct(t_mod *mod) |
---|
847 | { |
---|
848 | int i,j,k; |
---|
849 | triplet *triplet_struct; |
---|
850 | |
---|
851 | triplet_struct = (triplet *)mCalloc(1,sizeof(triplet)); |
---|
852 | triplet_struct->size = mod->ns; |
---|
853 | triplet_struct->pi_bc = (phydbl *)mCalloc(mod->ns,sizeof(phydbl )); |
---|
854 | triplet_struct->pi_cd = (phydbl *)mCalloc(mod->ns,sizeof(phydbl )); |
---|
855 | triplet_struct->pi_bd = (phydbl *)mCalloc(mod->ns,sizeof(phydbl )); |
---|
856 | triplet_struct->F_bc = (phydbl *)mCalloc(mod->ns*mod->ns*mod->ras->n_catg,sizeof(phydbl)); |
---|
857 | triplet_struct->F_cd = (phydbl *)mCalloc(mod->ns*mod->ns*mod->ras->n_catg,sizeof(phydbl)); |
---|
858 | triplet_struct->F_bd = (phydbl *)mCalloc(mod->ns*mod->ns,sizeof(phydbl)); |
---|
859 | triplet_struct->core = (phydbl ****)mCalloc(mod->ras->n_catg,sizeof(phydbl ***)); |
---|
860 | triplet_struct->p_one_site = (phydbl ***)mCalloc(mod->ns,sizeof(phydbl **)); |
---|
861 | triplet_struct->sum_p_one_site = (phydbl ***)mCalloc(mod->ns,sizeof(phydbl **)); |
---|
862 | triplet_struct->eigen_struct = (eigen *)Make_Eigen_Struct(mod->ns); |
---|
863 | triplet_struct->mod = mod; |
---|
864 | |
---|
865 | For(k,mod->ras->n_catg) |
---|
866 | { |
---|
867 | triplet_struct->core[k] = (phydbl ***)mCalloc(mod->ns,sizeof(phydbl **)); |
---|
868 | For(i,mod->ns) |
---|
869 | { |
---|
870 | triplet_struct->core[k][i] = (phydbl **)mCalloc(mod->ns,sizeof(phydbl *)); |
---|
871 | For(j,mod->ns) |
---|
872 | triplet_struct->core[k][i][j] = (phydbl *)mCalloc(mod->ns,sizeof(phydbl )); |
---|
873 | } |
---|
874 | } |
---|
875 | |
---|
876 | For(i,mod->ns) |
---|
877 | { |
---|
878 | triplet_struct->p_one_site[i] = (phydbl **)mCalloc(mod->ns,sizeof(phydbl *)); |
---|
879 | For(j,mod->ns) |
---|
880 | triplet_struct->p_one_site[i][j] = (phydbl *)mCalloc(mod->ns,sizeof(phydbl )); |
---|
881 | } |
---|
882 | |
---|
883 | For(i,mod->ns) |
---|
884 | { |
---|
885 | triplet_struct->sum_p_one_site[i] = (phydbl **)mCalloc(mod->ns,sizeof(phydbl *)); |
---|
886 | For(j,mod->ns) |
---|
887 | triplet_struct->sum_p_one_site[i][j] = (phydbl *)mCalloc(mod->ns,sizeof(phydbl )); |
---|
888 | } |
---|
889 | |
---|
890 | return triplet_struct; |
---|
891 | |
---|
892 | } |
---|
893 | |
---|
894 | ////////////////////////////////////////////////////////////// |
---|
895 | ////////////////////////////////////////////////////////////// |
---|
896 | |
---|
897 | void Make_Short_L(t_tree *tree) |
---|
898 | { |
---|
899 | if(!tree->short_l) |
---|
900 | tree->short_l = (phydbl *)mCalloc(tree->n_short_l,sizeof(phydbl)); |
---|
901 | } |
---|
902 | |
---|
903 | ////////////////////////////////////////////////////////////// |
---|
904 | ////////////////////////////////////////////////////////////// |
---|
905 | |
---|
906 | xml_attr *XML_Make_Attribute(xml_attr *prev, char *attr_name, char *attr_value) |
---|
907 | { |
---|
908 | xml_attr *new_attr; |
---|
909 | |
---|
910 | new_attr = (xml_attr *)mCalloc(1,sizeof(xml_attr)); |
---|
911 | |
---|
912 | new_attr->prev = prev; |
---|
913 | new_attr->next = NULL; |
---|
914 | if(prev != NULL) prev->next = new_attr; |
---|
915 | |
---|
916 | new_attr->name = (char *)mCalloc(strlen(attr_name)+1,sizeof(char)); |
---|
917 | strcpy(new_attr->name,attr_name); |
---|
918 | |
---|
919 | new_attr->value = (char *)mCalloc(strlen(attr_value)+1,sizeof(char)); |
---|
920 | strcpy(new_attr->value,attr_value); |
---|
921 | |
---|
922 | return new_attr; |
---|
923 | } |
---|
924 | |
---|
925 | ////////////////////////////////////////////////////////////// |
---|
926 | ////////////////////////////////////////////////////////////// |
---|
927 | |
---|
928 | void XML_Make_Node_Id(xml_node *n, char *id) |
---|
929 | { |
---|
930 | if(id) n->id = (char *)mCalloc(strlen(id)+1,sizeof(char)); |
---|
931 | } |
---|
932 | |
---|
933 | ////////////////////////////////////////////////////////////// |
---|
934 | ////////////////////////////////////////////////////////////// |
---|
935 | |
---|
936 | void XML_Make_Node_Value(xml_node *n, char *val) |
---|
937 | { |
---|
938 | if(val) n->value = (char *)mCalloc(strlen(val)+1,sizeof(char)); |
---|
939 | } |
---|
940 | |
---|
941 | ////////////////////////////////////////////////////////////// |
---|
942 | ////////////////////////////////////////////////////////////// |
---|
943 | |
---|
944 | |
---|
945 | xml_node *XML_Make_Node(char *name) |
---|
946 | { |
---|
947 | |
---|
948 | xml_node *new_node = (xml_node *)mCalloc(1,sizeof(xml_node)); |
---|
949 | |
---|
950 | if(name) |
---|
951 | { |
---|
952 | new_node->name = (char *)mCalloc(strlen(name)+1,sizeof(char)); |
---|
953 | } |
---|
954 | |
---|
955 | new_node->ds = (t_ds *)mCalloc(1,sizeof(t_ds)); |
---|
956 | |
---|
957 | return new_node; |
---|
958 | } |
---|
959 | |
---|
960 | ////////////////////////////////////////////////////////////// |
---|
961 | ////////////////////////////////////////////////////////////// |
---|
962 | |
---|
963 | void Make_Best_Spr(t_tree *tree) |
---|
964 | { |
---|
965 | tree->best_spr = Make_One_Spr(tree); |
---|
966 | Init_One_Spr(tree->best_spr); |
---|
967 | } |
---|
968 | |
---|
969 | ////////////////////////////////////////////////////////////// |
---|
970 | ////////////////////////////////////////////////////////////// |
---|
971 | |
---|
972 | void Make_Spr_List(t_tree *tree) |
---|
973 | { |
---|
974 | int i; |
---|
975 | |
---|
976 | tree->size_spr_list = 2*tree->n_otu-3; |
---|
977 | |
---|
978 | tree->spr_list = (t_spr **)mCalloc(2*tree->n_otu-2,sizeof(t_spr *)); |
---|
979 | |
---|
980 | For(i,2*tree->n_otu-2) |
---|
981 | { |
---|
982 | tree->spr_list[i] = Make_One_Spr(tree); |
---|
983 | Init_One_Spr(tree->spr_list[i]); |
---|
984 | } |
---|
985 | |
---|
986 | tree->perform_spr_right_away = NO; |
---|
987 | } |
---|
988 | |
---|
989 | ////////////////////////////////////////////////////////////// |
---|
990 | ////////////////////////////////////////////////////////////// |
---|
991 | |
---|
992 | t_spr *Make_One_Spr(t_tree *tree) |
---|
993 | { |
---|
994 | t_spr *a_spr; |
---|
995 | a_spr = (t_spr *)mCalloc(1,sizeof(t_spr)); |
---|
996 | a_spr->path = (t_node **)mCalloc(tree->n_otu,sizeof(t_node *)); |
---|
997 | return a_spr; |
---|
998 | } |
---|
999 | |
---|
1000 | ////////////////////////////////////////////////////////////// |
---|
1001 | ////////////////////////////////////////////////////////////// |
---|
1002 | |
---|
1003 | t_rate *RATES_Make_Rate_Struct(int n_otu) |
---|
1004 | { |
---|
1005 | t_rate *rates; |
---|
1006 | |
---|
1007 | rates = (t_rate *)mCalloc(1,sizeof(t_rate)); |
---|
1008 | rates->is_allocated = NO; |
---|
1009 | |
---|
1010 | if(n_otu > 0) |
---|
1011 | { |
---|
1012 | rates->is_allocated = YES; |
---|
1013 | rates->nd_r = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1014 | rates->br_r = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1015 | rates->buff_r = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1016 | rates->true_r = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1017 | rates->nd_t = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1018 | rates->buff_t = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1019 | rates->true_t = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1020 | rates->t_mean = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1021 | rates->t_prior = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1022 | rates->t_prior_min = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1023 | rates->t_prior_max = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1024 | rates->t_floor = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1025 | rates->t_ranked = (int *)mCalloc(2*n_otu-1,sizeof(int)); |
---|
1026 | rates->t_has_prior = (short int *)mCalloc(2*n_otu-1,sizeof(short int)); |
---|
1027 | rates->dens = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl)); |
---|
1028 | rates->triplet = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1029 | rates->n_jps = (int *)mCalloc(2*n_otu-1,sizeof(int)); |
---|
1030 | rates->t_jps = (int *)mCalloc(2*n_otu-2,sizeof(int)); |
---|
1031 | rates->cov_l = (phydbl *)mCalloc((2*n_otu-2)*(2*n_otu-2),sizeof(phydbl)); |
---|
1032 | rates->invcov = (phydbl *)mCalloc((2*n_otu-2)*(2*n_otu-2),sizeof(phydbl)); |
---|
1033 | rates->mean_l = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl)); |
---|
1034 | rates->ml_l = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl)); |
---|
1035 | rates->cur_l = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl)); |
---|
1036 | rates->u_ml_l = (phydbl *)mCalloc(2*n_otu-3,sizeof(phydbl)); |
---|
1037 | rates->u_cur_l = (phydbl *)mCalloc(2*n_otu-3,sizeof(phydbl)); |
---|
1038 | rates->cov_r = (phydbl *)mCalloc((2*n_otu-2)*(2*n_otu-2),sizeof(phydbl)); |
---|
1039 | rates->cond_var = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl)); |
---|
1040 | rates->mean_r = (phydbl *)mCalloc(2*n_otu-2,sizeof(phydbl)); |
---|
1041 | rates->mean_t = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1042 | rates->lca = (t_node **)mCalloc((2*n_otu-1)*(2*n_otu-1),sizeof(t_node *)); |
---|
1043 | rates->reg_coeff = (phydbl *)mCalloc((2*n_otu-2)*(2*n_otu-2),sizeof(phydbl)); |
---|
1044 | rates->trip_reg_coeff = (phydbl *)mCalloc((2*n_otu-2)*(6*n_otu-9),sizeof(phydbl)); |
---|
1045 | rates->trip_cond_cov = (phydbl *)mCalloc((2*n_otu-2)*9,sizeof(phydbl)); |
---|
1046 | rates->_2n_vect1 = (phydbl *)mCalloc(2*n_otu,sizeof(phydbl)); |
---|
1047 | rates->_2n_vect2 = (phydbl *)mCalloc(2*n_otu,sizeof(phydbl)); |
---|
1048 | rates->_2n_vect3 = (phydbl *)mCalloc(2*n_otu,sizeof(phydbl)); |
---|
1049 | rates->_2n_vect4 = (phydbl *)mCalloc(2*n_otu,sizeof(phydbl)); |
---|
1050 | rates->_2n_vect5 = (short int *)mCalloc(2*n_otu,sizeof(short int)); |
---|
1051 | rates->_2n2n_vect1 = (phydbl *)mCalloc(4*n_otu*n_otu,sizeof(phydbl)); |
---|
1052 | rates->_2n2n_vect2 = (phydbl *)mCalloc(4*n_otu*n_otu,sizeof(phydbl)); |
---|
1053 | rates->br_do_updt = (short int *)mCalloc(2*n_otu-1,sizeof(short int)); |
---|
1054 | rates->cur_gamma_prior_mean = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1055 | rates->cur_gamma_prior_var = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1056 | rates->n_tips_below = (int *)mCalloc(2*n_otu-1,sizeof(int)); |
---|
1057 | rates->time_slice_lims = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1058 | rates->n_time_slice_spans = (int *)mCalloc(2*n_otu-1,sizeof(int)); |
---|
1059 | rates->curr_slice = (int *)mCalloc(2*n_otu-1,sizeof(int)); |
---|
1060 | rates->has_survived = (int *)mCalloc(2*n_otu-1,sizeof(int)); |
---|
1061 | rates->survival_rank = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1062 | rates->survival_dur = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1063 | rates->calib_prob = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1064 | rates->curr_nd_for_cal = (int *)mCalloc(2*n_otu-1,sizeof(int)); |
---|
1065 | rates->t_prior_min_buff = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1066 | rates->t_prior_max_buff = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1067 | rates->times_partial_proba = (phydbl *)mCalloc(2*n_otu-1,sizeof(phydbl)); |
---|
1068 | } |
---|
1069 | |
---|
1070 | return rates; |
---|
1071 | } |
---|
1072 | |
---|
1073 | ////////////////////////////////////////////////////////////// |
---|
1074 | ////////////////////////////////////////////////////////////// |
---|
1075 | |
---|
1076 | t_cal *Make_Calib(int n_otu) |
---|
1077 | { |
---|
1078 | t_cal *calib; |
---|
1079 | int i; |
---|
1080 | i = 0; |
---|
1081 | calib = (t_cal *)mCalloc(1, sizeof(t_cal)); |
---|
1082 | calib -> proba = (phydbl *)mCalloc(2 * n_otu, sizeof(phydbl)); |
---|
1083 | calib -> all_applies_to = (t_node **)mCalloc(2 * n_otu - 1, sizeof(t_node *)); |
---|
1084 | For(i, 2 * n_otu - 1) |
---|
1085 | calib -> all_applies_to[i] = (t_node *)mCalloc(1, sizeof(t_node)); |
---|
1086 | return(calib); |
---|
1087 | } |
---|
1088 | |
---|
1089 | ////////////////////////////////////////////////////////////// |
---|
1090 | ////////////////////////////////////////////////////////////// |
---|
1091 | |
---|
1092 | void Make_Rmat_Weight(t_tree *mixt_tree) |
---|
1093 | { |
---|
1094 | t_tree *tree, *buff_tree; |
---|
1095 | |
---|
1096 | |
---|
1097 | tree = mixt_tree; |
---|
1098 | do |
---|
1099 | { |
---|
1100 | if(tree->is_mixt_tree == YES) tree = tree->next; |
---|
1101 | |
---|
1102 | buff_tree = mixt_tree->next; |
---|
1103 | do |
---|
1104 | { |
---|
1105 | if(buff_tree->mod->r_mat_weight == tree->mod->r_mat_weight) break; |
---|
1106 | buff_tree = buff_tree->next; |
---|
1107 | } |
---|
1108 | while(buff_tree != tree); |
---|
1109 | |
---|
1110 | if(buff_tree == tree) Free(tree->mod->r_mat_weight); |
---|
1111 | |
---|
1112 | tree = tree->next; |
---|
1113 | } |
---|
1114 | while(tree); |
---|
1115 | |
---|
1116 | |
---|
1117 | tree = mixt_tree; |
---|
1118 | do |
---|
1119 | { |
---|
1120 | if(tree->is_mixt_tree == YES) tree = tree->next; |
---|
1121 | tree->mod->r_mat_weight = NULL; |
---|
1122 | tree = tree->next; |
---|
1123 | } |
---|
1124 | while(tree); |
---|
1125 | |
---|
1126 | |
---|
1127 | tree = mixt_tree->next; |
---|
1128 | tree->mod->r_mat_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
1129 | Init_Scalar_Dbl(tree->mod->r_mat_weight); |
---|
1130 | tree->mod->r_mat_weight->v = 1.0; |
---|
1131 | |
---|
1132 | buff_tree = tree = mixt_tree; |
---|
1133 | do // For each mixt_tree |
---|
1134 | { |
---|
1135 | if(tree->is_mixt_tree == YES) tree = tree->next; |
---|
1136 | |
---|
1137 | buff_tree = mixt_tree->next; |
---|
1138 | do |
---|
1139 | { |
---|
1140 | if(buff_tree->mod->r_mat == tree->mod->r_mat) |
---|
1141 | { |
---|
1142 | tree->mod->r_mat_weight = buff_tree->mod->r_mat_weight; |
---|
1143 | break; |
---|
1144 | } |
---|
1145 | buff_tree = buff_tree->next; |
---|
1146 | } |
---|
1147 | while(buff_tree != tree); |
---|
1148 | |
---|
1149 | if(!tree->mod->r_mat_weight) |
---|
1150 | { |
---|
1151 | tree->mod->r_mat_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
1152 | Init_Scalar_Dbl(tree->mod->r_mat_weight); |
---|
1153 | tree->mod->r_mat_weight->v = 1.0; |
---|
1154 | } |
---|
1155 | |
---|
1156 | tree = tree->next; |
---|
1157 | |
---|
1158 | } |
---|
1159 | while(tree); |
---|
1160 | |
---|
1161 | tree = mixt_tree; |
---|
1162 | do |
---|
1163 | { |
---|
1164 | if(tree->next) |
---|
1165 | { |
---|
1166 | tree->mod->r_mat_weight->next = tree->next->mod->r_mat_weight; |
---|
1167 | tree->next->mod->r_mat_weight->prev = tree->mod->r_mat_weight; |
---|
1168 | } |
---|
1169 | else |
---|
1170 | { |
---|
1171 | tree->mod->r_mat_weight->next = NULL; |
---|
1172 | } |
---|
1173 | tree = tree->next; |
---|
1174 | } |
---|
1175 | while(tree); |
---|
1176 | |
---|
1177 | |
---|
1178 | } |
---|
1179 | |
---|
1180 | ////////////////////////////////////////////////////////////// |
---|
1181 | ////////////////////////////////////////////////////////////// |
---|
1182 | |
---|
1183 | void Make_Efrq_Weight(t_tree *mixt_tree) |
---|
1184 | { |
---|
1185 | t_tree *tree, *buff_tree; |
---|
1186 | |
---|
1187 | |
---|
1188 | tree = mixt_tree; |
---|
1189 | do |
---|
1190 | { |
---|
1191 | if(tree->is_mixt_tree == YES) tree = tree->next; |
---|
1192 | |
---|
1193 | buff_tree = mixt_tree->next; |
---|
1194 | do |
---|
1195 | { |
---|
1196 | if(buff_tree->mod->e_frq_weight == tree->mod->e_frq_weight) break; |
---|
1197 | buff_tree = buff_tree->next; |
---|
1198 | } |
---|
1199 | while(buff_tree != tree); |
---|
1200 | |
---|
1201 | if(buff_tree == tree) |
---|
1202 | { |
---|
1203 | Free(tree->mod->e_frq_weight); |
---|
1204 | } |
---|
1205 | |
---|
1206 | tree = tree->next; |
---|
1207 | } |
---|
1208 | while(tree); |
---|
1209 | |
---|
1210 | |
---|
1211 | tree = mixt_tree; |
---|
1212 | do |
---|
1213 | { |
---|
1214 | if(tree->is_mixt_tree == YES) tree = tree->next; |
---|
1215 | tree->mod->e_frq_weight = NULL; |
---|
1216 | tree = tree->next; |
---|
1217 | } |
---|
1218 | while(tree); |
---|
1219 | |
---|
1220 | |
---|
1221 | |
---|
1222 | tree = mixt_tree->next; |
---|
1223 | tree->mod->e_frq_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
1224 | Init_Scalar_Dbl(tree->mod->e_frq_weight); |
---|
1225 | tree->mod->e_frq_weight->v = 1.0; |
---|
1226 | |
---|
1227 | |
---|
1228 | buff_tree = tree = mixt_tree; |
---|
1229 | do // For each mixt_tree |
---|
1230 | { |
---|
1231 | if(tree->is_mixt_tree == YES) tree = tree->next; |
---|
1232 | |
---|
1233 | buff_tree = mixt_tree->next; |
---|
1234 | do |
---|
1235 | { |
---|
1236 | if(buff_tree->mod->e_frq == tree->mod->e_frq) |
---|
1237 | { |
---|
1238 | tree->mod->e_frq_weight = buff_tree->mod->e_frq_weight; |
---|
1239 | break; |
---|
1240 | } |
---|
1241 | buff_tree = buff_tree->next; |
---|
1242 | } |
---|
1243 | while(buff_tree != tree); |
---|
1244 | |
---|
1245 | if(!tree->mod->e_frq_weight) |
---|
1246 | { |
---|
1247 | tree->mod->e_frq_weight = (scalar_dbl *)mCalloc(1,sizeof(scalar_dbl)); |
---|
1248 | Init_Scalar_Dbl(tree->mod->e_frq_weight); |
---|
1249 | tree->mod->e_frq_weight->v = 1.0; |
---|
1250 | } |
---|
1251 | |
---|
1252 | tree = tree->next; |
---|
1253 | |
---|
1254 | } |
---|
1255 | while(tree); |
---|
1256 | |
---|
1257 | |
---|
1258 | tree = mixt_tree; |
---|
1259 | do |
---|
1260 | { |
---|
1261 | if(tree->next) |
---|
1262 | { |
---|
1263 | tree->mod->e_frq_weight->next = tree->next->mod->e_frq_weight; |
---|
1264 | tree->next->mod->e_frq_weight->prev = tree->mod->e_frq_weight; |
---|
1265 | } |
---|
1266 | else |
---|
1267 | { |
---|
1268 | tree->mod->e_frq_weight->next = NULL; |
---|
1269 | } |
---|
1270 | tree = tree->next; |
---|
1271 | } |
---|
1272 | while(tree); |
---|
1273 | } |
---|
1274 | |
---|
1275 | ////////////////////////////////////////////////////////////// |
---|
1276 | ////////////////////////////////////////////////////////////// |
---|
1277 | ////////////////////////////////////////////////////////////// |
---|
1278 | ////////////////////////////////////////////////////////////// |
---|
1279 | ////////////////////////////////////////////////////////////// |
---|
1280 | ////////////////////////////////////////////////////////////// |
---|