source: branches/lib/GDE/PHYML20130708/phyml/src/mpi_boot.c

Last change on this file was 17478, checked in by westram, 7 years ago
File size: 16.6 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 "mpi_boot.h"
14
15/* #ifdef MPI */
16
17/*********************************************************/
18
19void Bootstrap_MPI(t_tree *tree)
20{
21  int *site_num, n_site;
22  int replicate,j,k;
23  int position, init_len, nbRep;
24
25  calign *boot_data;
26  t_tree *boot_tree;
27  t_mod *boot_mod;
28  matrix *boot_mat;
29  char *s;
30
31  MPI_Status Stat;
32  int randomRecv, bootRecv, nbElem, i;
33  int *score_par, *score_tot;
34  char *bootStr, *t;
35 
36  randomRecv = nbElem = bootRecv = 0;
37
38  tree->print_boot_val       = 1;
39  tree->print_alrt_val       = 0;
40  boot_tree                  = NULL;
41
42  site_num = (int *)mCalloc(tree->data->init_len,sizeof(int));
43 
44  Free_Bip(tree);
45  Alloc_Bip(tree);
46  Get_Bip(tree->a_nodes[0],tree->a_nodes[0]->v[0],tree);
47
48  n_site = 0;
49  For(j,tree->data->crunch_len) For(k,tree->data->wght[j])
50    {
51      site_num[n_site] = j;
52      n_site++;
53    }
54
55  boot_data = Copy_Cseq(tree->data,tree->io);
56
57  if (Global_myRank == 0)
58    PhyML_Printf("\n. Non parametric bootstrap analysis \n");
59 
60  //number of bootstraps for each process
61  if (tree->mod->bootstrap%Global_numTask != 0) 
62    {
63      nbRep = (tree->mod->bootstrap / Global_numTask) + 1;
64      tree->mod->bootstrap = nbRep * Global_numTask;
65      if (Global_myRank == 0) {
66        PhyML_Printf("\n. The number of replicates is not a multiple of %d CPUs.\n", Global_numTask);
67        PhyML_Printf("\n. Will run %d replicates analysis.\n", tree->mod->bootstrap);
68      }
69    }
70  else
71    nbRep = tree->mod->bootstrap/Global_numTask;
72 
73  //Bip score
74  if (Global_myRank == 0) 
75    {
76      score_tot = (int *)mCalloc((2*tree->n_otu - 3),sizeof(int));
77      For(i,2*tree->n_otu-3)
78        score_tot[i] = 0;
79    }
80  else
81    score_tot = NULL;
82
83  score_par = (int *)mCalloc((2*tree->n_otu - 3),sizeof(int));
84  For(i,2*tree->n_otu-3)
85    score_par[i] = 0;
86
87  if (Global_myRank == 0)
88    PhyML_Printf("\n  [");
89
90  For(replicate, nbRep)
91    {
92      For(j,boot_data->crunch_len) boot_data->wght[j] = 0;
93
94      // Send random data to other process
95      if (Global_myRank == 0) {
96        // Compute number of data to send
97        if (tree->mod->bootstrap - randomRecv > Global_numTask)
98          nbElem = Global_numTask;
99        else
100          nbElem = tree->mod->bootstrap - randomRecv;
101
102        For(i,nbElem) {
103          For(j,boot_data->crunch_len) boot_data->wght[j] = 0;
104          init_len = 0;
105          // Create random data
106          For(j,boot_data->init_len)
107            {
108              position = Rand_Int(0,(int)(tree->data->init_len-1.0));
109              boot_data->wght[site_num[position]] += 1;
110              init_len++;
111            }
112           
113          Set_D_States(boot_data,tree->io->datatype,tree->io->state_len);
114
115          if (init_len != tree->data->init_len) {
116            MPI_Finalize();
117            Warn_And_Exit("\n== Pb. when copying sequences...\n");
118          }
119          // Send random data to other process, not to current process
120          if (i < nbElem-1) {
121            MPI_Ssend (boot_data->wght, boot_data->crunch_len, MPI_INT, i+1, Global_myRank, MPI_COMM_WORLD);
122#ifdef MPI_DEBUG
123fprintf (stderr, "\ntask %d, sending random to %d done\n", Global_myRank, i+1);
124fflush(stderr);
125#endif
126          }
127          randomRecv++;
128        }
129      }
130      else {
131        MPI_Recv (boot_data->wght, boot_data->crunch_len, MPI_INT, 0, MPI_ANY_TAG, MPI_COMM_WORLD, &Stat);
132#ifdef MPI_DEBUG
133fprintf (stderr, "\ntask %d, receiving random from task %d done\n", Global_myRank, Stat.MPI_SOURCE);
134fflush(stderr);
135#endif
136      }
137
138      init_len = 0;
139      For(j,boot_data->crunch_len) init_len += boot_data->wght[j];
140
141      if(init_len != tree->data->init_len) {
142        MPI_Finalize();
143        Warn_And_Exit("\n. Pb when copying sequences\n");
144      }
145
146      (tree->mod->io->datatype == NT)?
147        (Get_Base_Freqs(boot_data)):
148        (Get_AA_Freqs(boot_data));
149
150      if(tree->io->random_boot_seq_order) Randomize_Sequence_Order(boot_data);
151
152      boot_mod        = Copy_Model(tree->mod);
153      boot_mod->s_opt = tree->mod->s_opt; /* WARNING: re-using the same address here instead of creating a copying
154                                             requires to leave the value of s_opt unchanged during the bootstrap. */
155      boot_mod->io    = tree->io; /* WARNING: re-using the same address here instead of creating a copying
156                                     requires to leave the value of io unchanged during the bootstrap. */
157      Init_Model(boot_data,boot_mod,tree->io);
158
159      if(tree->io->mod->use_m4mod) M4_Init_Model(boot_mod->m4mod,boot_data,boot_mod);
160
161      if(tree->io->in_tree == 2)
162        {
163          switch(tree->io->tree_file_format)
164            {
165            case PHYLIP: 
166              {
167                rewind(tree->io->fp_in_tree);
168                boot_tree = Read_Tree_File_Phylip(tree->io->fp_in_tree);
169                break;
170              }
171            case NEXUS:
172              {
173                PhyML_Printf("\n. Unfortunately, PhyML cannot read NEXUS files and perform a bootstrap analysis."); 
174                PhyML_Printf("\n. Please use the PHYLIP format.."); 
175                PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
176                Warn_And_Exit("");
177                break;
178              }
179            default:
180              {
181                PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
182                Warn_And_Exit("");
183                break;
184              }
185            }
186        }
187      else
188        {
189          boot_mat = ML_Dist(boot_data,boot_mod);
190          boot_mat->tree = Make_Tree_From_Scratch(boot_data->n_otu,boot_data);
191          Fill_Missing_Dist(boot_mat);
192          Bionj(boot_mat);
193          boot_tree = boot_mat->tree;
194          boot_tree->mat = boot_mat;
195        }
196
197      boot_tree->mod                = boot_mod;
198      boot_tree->io                 = tree->io;
199      boot_tree->data               = boot_data;
200      boot_tree->mod->s_opt->print  = 0;
201      boot_tree->n_pattern          = boot_tree->data->crunch_len;
202      boot_tree->io->print_site_lnl = 0;
203      boot_tree->io->print_trace    = 0;
204
205      Set_Both_Sides(YES,boot_tree);
206     
207      if((boot_tree->mod->s_opt->random_input_tree) && (boot_tree->mod->s_opt->topo_search == SPR_MOVE)) Random_Tree(boot_tree);
208
209      Connect_CSeqs_To_Nodes(boot_data,boot_tree);
210      Share_Lk_Struct(tree,boot_tree);
211      Share_Spr_Struct(tree,boot_tree);
212      Share_Pars_Struct(tree,boot_tree);
213      Fill_Dir_Table(boot_tree);
214      Update_Dirs(boot_tree);
215
216      if(tree->mod->s_opt->greedy) Init_P_Lk_Tips_Double(boot_tree);
217      else                         Init_P_Lk_Tips_Int(boot_tree);
218      Init_Ui_Tips(boot_tree);
219      Init_P_Pars_Tips(boot_tree);
220      Br_Len_Not_Involving_Invar(boot_tree);
221     
222      if(boot_tree->mod->s_opt->opt_topo)
223        {
224          if(boot_tree->mod->s_opt->topo_search == NNI_MOVE) 
225            {
226              Simu_Loop(boot_tree);
227            }
228          else if((boot_tree->mod->s_opt->topo_search == SPR_MOVE) ||
229                  (boot_tree->mod->s_opt->topo_search == BEST_OF_NNI_AND_SPR))
230            {
231              Speed_Spr_Loop(boot_tree);
232            }
233        }
234      else
235        {
236          if(boot_tree->mod->s_opt->opt_subst_param || boot_tree->mod->s_opt->opt_bl)
237            Round_Optimize(boot_tree,boot_tree->data,ROUND_MAX);
238          else
239            Lk(NULL,boot_tree);
240        }
241
242      Free_Bip(boot_tree);
243
244      Alloc_Bip(boot_tree);
245
246      Match_Tip_Numbers(tree,boot_tree);
247
248      Get_Bip(boot_tree->a_nodes[0],
249              boot_tree->a_nodes[0]->v[0],
250              boot_tree);
251
252      Compare_Bip(tree,boot_tree,NO);
253     
254      Br_Len_Involving_Invar(boot_tree);
255
256      if(tree->io->print_boot_trees)
257        {
258          s = Write_Tree(boot_tree,NO);
259          t=(char *)mCalloc(T_MAX_LINE,sizeof(char));
260          Print_Fp_Out_Lines_MPI(boot_tree, tree->io, replicate+1, t);
261         
262          // Get bootstrap trees from other process and write to boot file
263          if (Global_myRank == 0) {
264            fprintf(tree->io->fp_out_boot_tree,"%s\n",s);
265            fprintf(tree->io->fp_out_boot_stats,"%s\n",t);
266            bootRecv++;
267            PhyML_Printf(".");
268            if(!((bootRecv)%tree->io->boot_prog_every)) 
269              {
270                PhyML_Printf("] %4d/%4d\n  ",bootRecv,tree->mod->bootstrap);
271                if(bootRecv != tree->mod->bootstrap) PhyML_Printf("[");
272              }
273           
274            // Compute number of bootstraps to receive
275            if (tree->mod->bootstrap - bootRecv > Global_numTask)
276              nbElem = Global_numTask;
277            else
278              nbElem = tree->mod->bootstrap - bootRecv + 1;
279             
280            bootStr=(char *)mCalloc(T_MAX_LINE,sizeof(char));
281            for (i=1; i<nbElem; i++) 
282              {
283                MPI_Recv (bootStr, T_MAX_LINE, MPI_CHAR, i, MPI_ANY_TAG, MPI_COMM_WORLD, &Stat);
284#ifdef MPI_DEBUG
285                PhyML_Fprintf (stderr, "\ntask %d, receiving bootstrap from task %d tag %d done\n", Global_myRank, Stat.MPI_SOURCE, Stat.MPI_TAG);
286                fflush(stderr);
287#endif
288                if (Stat.MPI_TAG == BootTreeTag)
289                  fprintf(tree->io->fp_out_boot_tree,"%s\n", bootStr);
290                if (Stat.MPI_TAG == BootStatTag)
291                  fprintf(tree->io->fp_out_boot_stats,"%s\n", bootStr);
292               
293                MPI_Recv (bootStr, T_MAX_LINE, MPI_CHAR, i, MPI_ANY_TAG, MPI_COMM_WORLD, &Stat);
294#ifdef MPI_DEBUG
295                fprintf (stderr, "\ntask %d, receiving bootstrap from task %d tag %d done\n", Global_myRank, Stat.MPI_SOURCE, Stat.MPI_TAG);
296                fflush(stderr);
297#endif
298                if (Stat.MPI_TAG == BootTreeTag)
299                  fprintf(tree->io->fp_out_boot_tree,"%s\n", bootStr);
300                if (Stat.MPI_TAG == BootStatTag)
301                  fprintf(tree->io->fp_out_boot_stats,"%s\n", bootStr);
302               
303                bootRecv++;
304                PhyML_Printf(".");
305                if(!((bootRecv)%tree->io->boot_prog_every)) {
306                  PhyML_Printf("] %4d/%4d\n  ",bootRecv,tree->mod->bootstrap);
307                  if(bootRecv != tree->mod->bootstrap) PhyML_Printf("[");
308                }
309              }
310            Free(bootStr);
311          }
312          else {
313            MPI_Ssend (s, T_MAX_LINE, MPI_CHAR, 0, BootTreeTag, MPI_COMM_WORLD);
314            MPI_Ssend (t, T_MAX_LINE, MPI_CHAR, 0, BootStatTag, MPI_COMM_WORLD);
315#ifdef MPI_DEBUG
316fprintf (stderr, "\ntask %d, sending bootstraps done\n", Global_myRank);
317fflush(stderr);
318#endif
319          }
320          Free (t);
321          Free(s);
322        }
323
324#ifndef QUIET
325fflush(stdout);
326#endif
327
328      if(boot_tree->mat) Free_Mat(boot_tree->mat);
329      Free_Tree(boot_tree);     
330      Free_Model(boot_mod);
331
332      //Each process computes the Bip score sum for all its bootstrap trees
333      For(i,2*tree->n_otu-3) score_par[i] = tree->a_edges[i]->bip_score;
334
335      //Each process sends its Bip score sum. The sums are summed.
336      MPI_Reduce(score_par, score_tot, 2*tree->n_otu - 3, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
337    }
338
339 
340  if (Global_myRank == 0) 
341    {
342      For(i,2*tree->n_otu-3)
343        tree->a_edges[i]->bip_score = score_tot[i];
344      Free (score_tot);
345    }
346  Free (score_par);
347
348  if (Global_myRank == 0)
349    if(((bootRecv)%tree->io->boot_prog_every)) PhyML_Printf("] %4d/%4d\n ",bootRecv,tree->mod->bootstrap);
350
351  tree->lock_topo = 1; /* TopoLOGy should not be modified afterwards */
352
353  if(tree->io->print_boot_trees)
354    {
355      fclose(tree->io->fp_out_boot_tree);
356      fclose(tree->io->fp_out_boot_stats);
357    }
358
359  Free_Cseq(boot_data);
360  Free(site_num);
361}
362
363/*********************************************************/
364
365void Print_Fp_Out_Lines_MPI(t_tree *tree, option *io, int n_data_set, char *bootStr)
366{
367  char *s, *tmp;
368
369  // Build a string to be sent to the writing process
370  s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
371  tmp=(char *)mCalloc(T_MAX_LINE,sizeof(char));
372 
373  if (Global_myRank == 0 && n_data_set == 1) {
374    snprintf(tmp, T_MAX_LINE, ". Sequence file : [%s]\n\n", Basename(io->in_align_file)); strncat (s, tmp, T_MAX_LINE);
375   
376    (tree->mod->io->datatype == NT)?
377      (snprintf(tmp, T_MAX_LINE, ". Model of nucleotides substitution : %s\n\n", io->mod->modelname->s)):
378      (snprintf(tmp, T_MAX_LINE, ". Model of amino acids substitution : %s\n\n", io->mod->modelname->s));
379    strncat (s, tmp, T_MAX_LINE);
380   
381    switch(io->in_tree)
382      {
383      case 0: { snprintf(tmp, T_MAX_LINE, ". Initial tree : [BioNJ]\n\n");               break; }
384      case 1: { snprintf(tmp, T_MAX_LINE, ". Initial tree : [parsimony]\n\n");           break; }
385      case 2: { snprintf(tmp, T_MAX_LINE, ". Initial tree : [%s]\n\n",io->in_tree_file); break; }
386      }
387    strncat (s, tmp, T_MAX_LINE);
388
389    strncat (s, "\n", T_MAX_LINE);
390   
391    /*headline 1*/
392    strncat (s, ". Data\t", T_MAX_LINE);
393
394    strncat (s, "Nb of \t", T_MAX_LINE);
395
396    strncat (s, "Likelihood\t", T_MAX_LINE);
397
398    strncat (s, "Discrete   \t", T_MAX_LINE);
399
400    if(tree->mod->ras->n_catg > 1)
401      strncat (s, "Number of \tGamma shape\t", T_MAX_LINE);
402
403    strncat (s, "Proportion of\t", T_MAX_LINE);
404
405    if(tree->mod->whichmodel <= 6)
406      strncat (s, "Transition/ \t", T_MAX_LINE);
407   
408    strncat (s, "Nucleotides frequencies               \t", T_MAX_LINE);
409
410    if((tree->mod->whichmodel == GTR) ||
411       (tree->mod->whichmodel == CUSTOM))
412      strncat (s, "Instantaneous rate matrix              \t", T_MAX_LINE);
413
414    strncat (s, "\n", T_MAX_LINE);
415
416    /*headline 2*/
417    strncat (s, "  set\t", T_MAX_LINE);
418
419    strncat (s, "taxa\t", T_MAX_LINE);
420
421    strncat (s, "LOGlk     \t", T_MAX_LINE);
422
423    strncat (s, "gamma model\t", T_MAX_LINE);
424
425    if(tree->mod->ras->n_catg > 1)
426      strncat (s, "categories\tparameter  \t", T_MAX_LINE);
427
428    strncat (s, "invariant    \t", T_MAX_LINE);
429
430    if(tree->mod->whichmodel <= 6)
431      strncat (s, "transversion\t", T_MAX_LINE);
432
433    strncat (s, "f(A)      f(C)      f(G)      f(T)    \t", T_MAX_LINE);
434
435    if((tree->mod->whichmodel == GTR) ||
436       (tree->mod->whichmodel == CUSTOM))
437      strncat (s, "[A---------C---------G---------T------]\t", T_MAX_LINE);
438
439    strncat (s, "\n", T_MAX_LINE);
440
441    /*headline 3*/
442    if(tree->mod->whichmodel == TN93) {
443      strncat (s, "    \t      \t          \t           \t", T_MAX_LINE);
444      if(tree->mod->ras->n_catg > 1)
445        strncat (s, "         \t         \t", T_MAX_LINE);
446      strncat (s, "             \t", T_MAX_LINE);
447      strncat (s, "purines pyrimid.\t", T_MAX_LINE);
448      strncat (s, "\n", T_MAX_LINE);
449    }
450
451    strncat (s, "\n", T_MAX_LINE);
452  }
453
454  /*line items*/
455
456  snprintf(tmp, T_MAX_LINE, "  #%d\t", (((n_data_set-1)*Global_numTask)+Global_myRank+1));
457  strncat (s, tmp, T_MAX_LINE);
458 
459  snprintf(tmp, T_MAX_LINE, "%d   \t",tree->n_otu); strncat (s, tmp, T_MAX_LINE);
460 
461  snprintf(tmp, T_MAX_LINE, "%.5f\t",tree->c_lnL); strncat (s, tmp, T_MAX_LINE);
462 
463  snprintf(tmp, T_MAX_LINE, "%s        \t",
464          (tree->mod->ras->n_catg>1)?("Yes"):("No ")); strncat (s, tmp, T_MAX_LINE);
465 
466  if(tree->mod->ras->n_catg > 1)
467    {
468      snprintf(tmp, T_MAX_LINE, "%d        \t",tree->mod->ras->n_catg); strncat (s, tmp, T_MAX_LINE);
469      snprintf(tmp, T_MAX_LINE, "%.3f    \t",tree->mod->ras->alpha->v); strncat (s, tmp, T_MAX_LINE);
470    }
471 
472  snprintf(tmp, T_MAX_LINE, "%.3f    \t",tree->mod->ras->pinvar->v); strncat (s, tmp, T_MAX_LINE);
473 
474  if(tree->mod->whichmodel <= 5)
475    {
476      snprintf(tmp, T_MAX_LINE, "%.3f     \t",tree->mod->kappa->v); strncat (s, tmp, T_MAX_LINE);
477    }
478  else if(tree->mod->whichmodel == TN93)
479    {
480      snprintf(tmp, T_MAX_LINE, "%.3f   ",
481              tree->mod->kappa->v*2.*tree->mod->lambda->v/(1.+tree->mod->lambda->v));
482      strncat (s, tmp, T_MAX_LINE);
483      snprintf(tmp, T_MAX_LINE, "%.3f\t",
484              tree->mod->kappa->v*2./(1.+tree->mod->lambda->v));
485      strncat (s, tmp, T_MAX_LINE);
486    }
487
488  if(tree->mod->io->datatype == NT)
489    {
490      snprintf(tmp, T_MAX_LINE, "%8.5f  ",tree->mod->e_frq->pi->v[0]); strncat (s, tmp, T_MAX_LINE);
491      snprintf(tmp, T_MAX_LINE, "%8.5f  ",tree->mod->e_frq->pi->v[1]); strncat (s, tmp, T_MAX_LINE);
492      snprintf(tmp, T_MAX_LINE, "%8.5f  ",tree->mod->e_frq->pi->v[2]); strncat (s, tmp, T_MAX_LINE);
493      snprintf(tmp, T_MAX_LINE, "%8.5f\t",tree->mod->e_frq->pi->v[3]); strncat (s, tmp, T_MAX_LINE);
494    }
495   
496  if((tree->mod->whichmodel == GTR) || (tree->mod->whichmodel == CUSTOM))
497    {
498      int i,j;
499
500      For(i,4)
501        {
502          if (i!=0) {
503            /*format*/
504            snprintf(tmp, T_MAX_LINE, "      \t     \t          \t           \t");
505            strncat (s, tmp, T_MAX_LINE);
506            if(tree->mod->ras->n_catg > 1) {
507              snprintf(tmp, T_MAX_LINE, "          \t           \t");
508              strncat (s, tmp, T_MAX_LINE);
509            }
510            snprintf(tmp, T_MAX_LINE, "             \t                                      \t");
511            strncat (s, tmp, T_MAX_LINE);
512          }
513          For(j,4) {
514            snprintf(tmp, T_MAX_LINE, "%8.5f  ",tree->mod->r_mat->qmat->v[i*4+j]);
515            strncat (s, tmp, T_MAX_LINE);
516          }
517          if (i<3) {
518            snprintf(tmp, T_MAX_LINE, "\n");
519            strncat (s, tmp, T_MAX_LINE);
520          }
521        }
522    }
523   
524    Free (tmp);
525    strncpy (bootStr, s, T_MAX_LINE);
526    Free (s);
527
528  return;
529}
530/* #endif */
Note: See TracBrowser for help on using the repository browser.