source: branches/ali/GDE/PHYML20130708/phyml/src/make.c

Last change on this file was 10307, checked in by aboeckma, 11 years ago

added most recent version of phyml

File size: 39.8 KB
Line 
1/*
2
3PhyML:  a program that  computes maximum likelihood phylogenies from
4DNA or AA homologous sequences.
5
6Copyright (C) Stephane Guindon. Oct 2003 onward.
7
8All parts of the source except where indicated are distributed under
9the GNU public licence. See http://www.opensource.org for details.
10
11*/
12
13#include "make.h"
14
15//////////////////////////////////////////////////////////////
16
17void 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
52void 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
66void 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
85void 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
105t_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->= 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
152void 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
161void 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
171void 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
181void 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
200void 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
247void 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
293void 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
305nni *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
316t_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
338void Make_Node_Lk(t_node *n)
339{
340/*   n->n_ex_nodes = (int *)mCalloc(2,sizeof(int)); */
341  return;
342}
343
344nexcom **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
365nexparm *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
379matrix *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
410t_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
431t_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
443void 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
452void 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
469void 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
479calign *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
515t_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
529t_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
539void 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
563t_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
576t_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
622void 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
672t_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 */
710void 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
725t_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
745t_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
779option *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
810t_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
823eigen *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
846triplet *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
897void 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
906xml_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
928void 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
936void 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
945xml_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
963void 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
972void 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
992t_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
1003t_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
1076t_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
1092void 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
1183void 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//////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.