source: trunk/GDE/PHYML/main.c

Last change on this file was 19480, checked in by westram, 17 months ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.5 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 "utilities.h"
14#include "lk.h"
15#include "optimiz.h"
16#include "bionj.h"
17#include "models.h"
18#include "free.h"
19#include "options.h"
20#include "simu.h"
21
22
23#ifdef PHYML
24
25/* int T_MAX_FILE; */
26
27int main(int argc, char **argv)
28{
29  seq **data;
30  allseq *alldata;
31  option *input;
32  char *s_tree, *s_any;
33  FILE *fp_phyml_tree,*fp_phyml_stats,*fp_phyml_lk,*fp_best_tree,*fp_best_tree_stats;
34  arbre *tree;
35  int n_data_sets;
36  matrix *mat;
37  model *mod;
38  time_t t_beg,t_end;
39  div_t hour,min;
40  int num_tree,tree_line_number;
41  double best_loglk;
42 
43
44  srand(time(NULL));
45
46  tree = NULL;
47  mod  = NULL;
48
49/*   Init_Constant(); */
50
51  s_any = (char *)mCalloc(T_MAX_FILE,sizeof(char));
52
53  fflush(NULL);
54
55  input = (option *)Get_Input(argc,argv);
56
57  Make_Model_Complete(input->mod);
58
59  mod = input->mod;
60
61  fp_phyml_stats = Openfile(input->phyml_stat_file,input->phyml_stat_file_open_mode);
62
63  fprintf(fp_phyml_stats,"\n- PHYML %s -\n\n", VERSION);
64
65  fp_phyml_tree = Openfile(input->phyml_tree_file,input->phyml_tree_file_open_mode);
66
67  n_data_sets = 0;
68
69  fp_phyml_lk = fopen(input->phyml_lk_file,"w");
70
71  if(input->inputtree) Test_Multiple_Data_Set_Format(input);
72  else input->n_trees = 1;
73 
74/*   if(input->n_data_sets > 1) input->n_trees = 1; */
75
76  best_loglk = UNLIKELY;
77  tree_line_number = 0;
78
79  do
80    {
81
82      n_data_sets++;
83
84      time(&t_beg);
85
86      if(n_data_sets > input->n_data_sets) 
87        {
88          data = NULL;   
89        }
90      else
91        {
92            data = Get_Seq(input);
93        }
94
95      if(data)
96        {
97          if(n_data_sets > 1) printf("\n. Data set [#%d]\n",n_data_sets);
98
99          printf("\n. Compressing sequences...\n");
100
101          alldata = Compact_Seq(data,input);
102
103          Free_Seq(data,alldata->n_otu);
104
105          Init_Model(alldata,mod);
106         
107          Check_Ambiguities(alldata,input->mod->datatype,input->mod->stepsize);
108
109
110          For(num_tree,input->n_trees)
111            {
112              if(!input->inputtree)
113                {
114                  printf("\n. Computing pairwise distances...\n");
115                 
116                  mat = ML_Dist(alldata,mod);
117
118                  printf("\n. Building BIONJ tree...\n");
119                             
120                  mat->tree = Make_Tree(alldata);
121
122                  Bionj(mat);
123                                 
124                  tree = mat->tree;
125
126                  Free_Mat(mat);
127                }
128              else 
129                {
130                    if(input->n_trees > 1) printf("\n. Reading user tree [#%d]\n",tree_line_number+1);
131                    else printf("\n. Reading user tree...\n");
132                 
133                    if(input->n_trees == 1) 
134                        {
135                            rewind(input->fp_input_tree);
136                            tree_line_number = 0;
137                        }
138
139                    tree = Read_Tree_File(input->fp_input_tree);
140                   
141                    tree_line_number++;
142                 
143                    if(!tree) 
144                        {
145                            printf("\n. Missing tree for data set #%d\n",n_data_sets);
146                            printf("  This data set is not analyzed.\n");
147                            data = NULL;
148                        }
149                   
150                    if(!tree->has_branch_lengths)
151                        {
152                            printf("\n. Computing branch length estimates...\n");
153                           
154                            Order_Tree_CSeq(tree,alldata);
155                           
156                            mat = ML_Dist(alldata,mod);
157                           
158                            mat->tree = tree;
159                     
160                            mat->method = 0;
161                           
162                            Bionj_Br_Length(mat);
163                           
164                            Free_Mat(mat);
165                        }
166                }
167             
168              if(!tree) continue;
169             
170             
171              tree->mod        = mod;
172              tree->input      = input;
173              tree->data       = alldata;
174              tree->both_sides = 1;
175              tree->n_pattern  = tree->data->crunch_len/tree->mod->stepsize;
176             
177              Order_Tree_CSeq(tree,alldata);
178             
179              Make_Tree_4_Lk(tree,alldata);
180             
181              if(tree->mod->s_opt->opt_topo)
182                Simu(tree,1000);
183              else
184                {
185                  if(tree->mod->s_opt->opt_free_param)
186                    Round_Optimize(tree,tree->data);
187                  else
188                    {
189                      Lk(tree,tree->data);
190                      printf("\n. Log(lk) :               * -> %15.6f ",tree->tot_loglk);
191                    }
192                }
193             
194              if(tree->mod->bootstrap) Bootstrap(tree);
195             
196              Update_BrLen_Invar(tree);
197             
198              s_tree = Write_Tree(tree);
199             
200              fprintf(fp_phyml_tree,"%s\n",s_tree);
201
202              Free(s_tree);
203
204              Unconstraint_Lk(tree);
205             
206              time(&t_end);
207             
208              hour = div(t_end-t_beg,3600);
209              min  = div(t_end-t_beg,60  );
210             
211              min.quot -= hour.quot*60;
212             
213              if (input->n_data_sets==1)
214                Print_Fp_Out(fp_phyml_stats, t_beg, t_end, tree, input, n_data_sets);
215              else
216                Print_Fp_Out_Lines(fp_phyml_stats, t_beg, t_end, tree, input, n_data_sets);
217             
218              printf("\n\n. Time used %dh%dm%ds\n", hour.quot,min.quot,(int)(t_end-t_beg)%60);
219              printf("\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n");
220             
221              Print_Site_Lk(tree,fp_phyml_lk);
222/*            fprintf(fp_phyml_lk,"%f\n",tree->tot_loglk);        */
223             
224             
225              if((input->n_data_sets == 1) &&
226                 (input->n_trees     >  1) &&
227                 (tree->tot_loglk > best_loglk))
228                {
229
230                  best_loglk = tree->tot_loglk;
231                  strcpy(s_any,input->seqfile);
232                  fp_best_tree = fopen(s_any = strcat(s_any,"_phyml_best_tree.txt"),"w");
233                  s_tree = Write_Tree(tree);
234                  fprintf(fp_best_tree,"%s\n",s_tree);
235                  Free(s_tree);
236                  strcpy(s_any,input->seqfile);
237                  fp_best_tree_stats = fopen(s_any = strcat(s_any,"_phyml_best_stat.txt"),"w");
238                  Print_Fp_Out(fp_best_tree_stats, 
239                               t_beg, 
240                               t_end, 
241                               tree, 
242                               input, 
243                               n_data_sets);
244                 
245                  fclose(fp_best_tree);
246                  fclose(fp_best_tree_stats);
247                }
248
249              Free_Tree_Lk(tree);
250             
251              Free_Tree(tree);
252             
253              if(input->n_data_sets > 1) 
254                  {
255                      break;
256                  }
257            }
258          Free_Cseq(alldata);
259        }
260    }while(data);
261
262  Free_Model(mod);
263
264  if(input->fp_seq ) fclose(input->fp_seq );
265  if(input->fp_input_tree) fclose(input->fp_input_tree);
266
267  fclose(fp_phyml_lk);
268
269  fclose(fp_phyml_tree);
270
271  fclose(fp_phyml_stats);
272
273  Free_Input(input);   
274 
275  Free(s_any);
276
277  return 0;
278}
279
280#endif
281/*********************************************************/
282
Note: See TracBrowser for help on using the repository browser.