~ubuntu-branches/ubuntu/raring/paml/raring

« back to all changes in this revision

Viewing changes to Technical/Simulation/multiruns.c

  • Committer: Bazaar Package Importer
  • Author(s): Pjotr Prins
  • Date: 2010-09-11 23:01:37 UTC
  • Revision ID: james.westby@ubuntu.com-20100911230137-jjf5d0blx5p0m9ba
Tags: upstream-4.4c
ImportĀ upstreamĀ versionĀ 4.4c

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* multiruns.c
 
2
 
 
3
     cl -O2 -Ot -W4 multiruns.c
 
4
     cc -o multiruns -O2 multiruns.c -lm
 
5
 
 
6
     multiruns <rstfile1> <rstfile2> ... <lnLColumn> 
 
7
 
 
8
  Examples (comparing three runs with lnL in column 19 in rst1):
 
9
 
 
10
     multiruns rst1.r1 rst1.r2 rst1.r3 19
 
11
     multiruns a/rst1 b/rst1 c/rst1 19
 
12
 
 
13
   March 2003, Ziheng Yang
 
14
   September 2005, changed tworuns into multiruns, ziheng yang
 
15
 
 
16
   This program compares outputs from multiple separate ML runs analyzing
 
17
   many data sets (using ndata) to assemble a result file.  Because of local 
 
18
   peaks and convergence problems, multiple runs for the same analysis may not 
 
19
   generate the same results.  Then we should use the results corresponding to 
 
20
   the highest lnL.  This program takes input files which have summary results 
 
21
   from multiple runs, one line for each data set.  The program takes one line 
 
22
   from each of the input files and compare the first field, which is an index 
 
23
   column and should be identical between the input files, and an lnL column.  
 
24
   The program decides which run generated the highest lnL, and copy the line 
 
25
   from that run into the output file: out.txt.
 
26
 
 
27
   This is useful when you analyze the same set of simulated replicate data 
 
28
   sets multiple times, using different starting values.  For example, codeml 
 
29
   may write a line of output in rst1 for each data set, including parameter 
 
30
   estimates and lnL.  You can then use this program to compare the rst1 output 
 
31
   files from multiple runs to generate one output file.  The program allows the 
 
32
   fields to be either numerical or text, but the first (index) and lnL columns
 
33
   should be numerical.
 
34
 
 
35
*/
 
36
 
 
37
#include <stdio.h>
 
38
#include <stdlib.h>
 
39
#include <string.h>
 
40
#include <ctype.h>
 
41
 
 
42
#define MAXNFIILES  20
 
43
#define MAXNFIELDS  1000
 
44
#define MAXLLINE    64000
 
45
 
 
46
void error2 (char * message);
 
47
FILE *gfopen(char *filename, char *mode);
 
48
int splitline (char line[], int fields[]);
 
49
 
 
50
 
 
51
int main(int argc, char* argv[])
 
52
{
 
53
   FILE *fout, *fin[MAXNFIILES];
 
54
   char infile[MAXNFIILES][96]={"rst1.r1", "rst1.r2"}, outfile[96]="out.txt";
 
55
   int index=0, nfile, nfileread, lnLcolumn=13, i, nrecords=0, lline=MAXLLINE;
 
56
   int nfields[MAXNFIILES],fields[MAXNFIELDS], minf, maxf, miss[MAXNFIILES];
 
57
   char *line[MAXNFIILES];
 
58
   double lnL[MAXNFIILES], lnLmin, lnLmax, indexfield[MAXNFIILES], y;
 
59
 
 
60
   puts("Usage: \n\tmultiruns <file1> <file2> ... <lnLcolumn>\n");
 
61
   nfile=argc-2;
 
62
   if(nfile<2) exit(1);
 
63
 
 
64
   for(i=0; i<nfile; i++) { 
 
65
      strcpy(infile[i], argv[1+i]);
 
66
      fin[i]=gfopen(infile[i],"r");
 
67
      if((line[i]=(char*)malloc(MAXLLINE*sizeof(char)))==NULL) error2("oom");
 
68
      printf("%s  ", infile[i]);
 
69
   }
 
70
   printf("  ==>  %s\n\n", outfile);
 
71
 
 
72
   sscanf(argv[argc-1], "%d", &lnLcolumn);
 
73
   lnLcolumn--;
 
74
   fout=(FILE*)gfopen(outfile,"w");
 
75
 
 
76
   for(nrecords=0; ; nrecords++) {
 
77
      for(i=0,nfileread=0; i<nfile; i++) { 
 
78
         nfields[i]=0; lnL[i]=0; line[i][0]='\0';  miss[i]=1;
 
79
         if(!fgets(line[i], lline, fin[i])) continue;
 
80
         nfields[i]=splitline(line[i],fields);
 
81
 
 
82
         if(nfields[i]>index) sscanf(line[i]+fields[index], "%lf", &indexfield[i]);
 
83
         if(nfields[i]>lnLcolumn) {
 
84
            sscanf(line[i]+fields[lnLcolumn], "%lf", &lnL[i]);
 
85
            miss[i]=0;
 
86
            nfileread++;
 
87
         }
 
88
      }
 
89
      if(nfileread==0) break;
 
90
      for(i=0,y=-1; i<nfile; i++) {
 
91
         if(!miss[i]) {
 
92
           if(y==-1) y=indexfield[i];
 
93
           else if(y!=indexfield[i]) error2("index field different");
 
94
         }
 
95
      }
 
96
      for(i=0,lnLmin=1e300,lnLmax=-1e300; i<nfile; i++) {
 
97
         if(miss[i]) continue;
 
98
         if(lnL[i]<lnLmin) { lnLmin=lnL[i];  minf=i; }
 
99
         if(lnL[i]>lnLmax) { lnLmax=lnL[i];  maxf=i; }
 
100
      }
 
101
 
 
102
      fprintf(fout, "%s", line[maxf]);
 
103
 
 
104
      printf("record %4d  (", nrecords+1);
 
105
      for(i=0; i<nfile; i++)  printf("%c", (miss[i]?'-':'+')); 
 
106
      printf(") %10.3f (%d) - %10.3f (%d) = %8.3f\n", 
 
107
         lnLmin, minf+1, lnLmax, maxf+1, lnLmin-lnLmax);
 
108
   }
 
109
   printf("\nwrote %d records into %s\n", nrecords, outfile);
 
110
   for(i=0; i<nfile; i++) { fclose(fin[i]);  free(line[i]); }
 
111
   fclose(fout);
 
112
}
 
113
 
 
114
void error2 (char * message)
 
115
{ printf("\nError: %s.\n", message); exit(-1); }
 
116
 
 
117
FILE *gfopen(char *filename, char *mode)
 
118
{
 
119
   FILE *fp;
 
120
 
 
121
   if(filename==NULL || filename[0]==0) 
 
122
      error2("file name empty.");
 
123
 
 
124
   fp=(FILE*)fopen(filename, mode);
 
125
   if(fp==NULL) {
 
126
      printf("\nerror when opening file %s\n", filename);
 
127
      if(!strchr(mode,'r')) exit(-1);
 
128
      printf("tell me the full path-name of the file? ");
 
129
      scanf("%s", filename);
 
130
      if((fp=(FILE*)fopen(filename, mode))!=NULL)  return(fp);
 
131
      puts("Can't find the file.  I give up.");
 
132
      exit(-1);
 
133
   }
 
134
   return(fp);
 
135
}
 
136
 
 
137
int splitline (char line[], int fields[])
 
138
{
 
139
/* This finds out how many fields there are in the line, and marks the starting 
 
140
   positions of the fields.
 
141
   Fields are separated by spaces, and texts are allowed as well.
 
142
*/
 
143
   int lline=64000, i, nfields=0, InSpace=1;
 
144
   char *p=line;
 
145
 
 
146
   for(i=0; i<lline && *p && *p!='\n'; i++,p++) {
 
147
      if (isspace(*p))
 
148
         InSpace=1;
 
149
      else  {
 
150
         if(InSpace) {
 
151
            InSpace=0;
 
152
            fields[nfields++]=i;
 
153
            /* if(nfields>MAXNFIELDS) puts("raise MAXNFIELDS?"); */
 
154
         }
 
155
      }
 
156
   }
 
157
   return(nfields);
 
158
}
 
159