1 /*------------------------------------------------------
2 Maximum likelihood estimation
3 of migration rate and effectice population size
4 using a Metropolis-Hastings Monte Carlo algorithm
5 -------------------------------------------------------
6 W O R L D R O U T I N E S
7
8 creates tree structures,
9 reads tree [has to be done
10
11 prints results,
12 and finally helps to destroy itself.
13
14 Peter Beerli 1996, Seattle
15 beerli@genetics.washington.edu
16 $Id$
17
18 modifed for use with recombine by Jon Yamato (1/8/97)
19 -------------------------------------------------------*/
20
21 #include "world.h"
22
23 #ifndef REC_MODELLIKE_INCLUDE
24 #include "rec_modellike.h"
25 #endif
26
27 #ifdef DMALLOC_FUNC_CHECK
28 #include "dmalloc.h"
29 #endif
30
31 #ifdef DMEMDEBUG
32 #define MEMDEBUG
33 #include "memdebug.h"
34 #endif
35
36 #define PLANESIZE 36
37 #define PLANEBIGTICKS 6
38 #define PLANEBIGTICKVALUES {-3, -2, -1, 0, 1, 2}
39 #define PLANETICKS {0.001, 0.0013895, 0.0019307, 0.0026827, 0.00372759, \
40 0.00517947, 0.00719686, 0.01, 0.013895, 0.019307, \
41 0.026827, 0.0372759, 0.0517947, 0.0719686, 0.1, \
42 0.13895, 0.19307, 0.26827, 0.372759, 0.517947, \
43 0.719686, 1., 1.3895, 1.9307, 2.6827, 3.72759, \
44 5.17947, 7.19686, 10., 13.895, 19.307, 26.827, 37.2759, \
45 51.7947, 71.9686, 100.}
46 #define CONTOURLEVELS 8
47 #define CONTOURS_LOCUS {0.0,-1.38629/2., -5.99147/2., -9.21034/2.,\
48 0.0,-1.38629/2., -5.99147/2., -9.21034/2.}
49 #if 0
50 #define CONTOURS_LOCUS {0.0,-3.35670/2., -9.48773/2., -13.2767/2.,\
51 0.0,-3.35670/2., -9.48773/2., -13.2767/2.}
52 #define CONTOURS_LOCI {0.0,-4.35146/2., -11.0705/2., -15.0863/2.,\
53 0.0,-4.35146/2., -11.0705/2., -15.0863/2.}
54 #endif
55
56 #define PAGEFEED fprintf(outfile,"\n\n\n\n");
57
58 extern long totchains;
59 extern FILE *outfile;
60
61 extern double model_likelihood(option_struct *op, data_fmt *data, double theta,
62 double rec, long lowcus, long chain);
63 extern double fmodel_likelihood(option_struct *op, data_fmt *data, double theta,
64 double rec, double **denom, long lowcus, long chain);
65
66
plot_surface(option_struct * op,data_fmt * data,long lowcus,char *** plane,long x)67 void plot_surface(option_struct *op, data_fmt *data, long lowcus, char ***plane, long x)
68 {
69 long i;
70 long loci;
71
72 loci = getdata_numloci(op,data);
73
74 fprintf(outfile, "\nLog-Likelihood surface\n");
75 fprintf(outfile, "---------------------------------------\n\n");
76 fprintf(outfile, "Legend:\n");
77 fprintf(outfile, " X = Maximum likelihood\n");
78 fprintf(outfile, " * = in approximative 50%% confidence limit\n");
79 fprintf(outfile, " + = in approximative 95%% confidence limit\n");
80 fprintf(outfile, " - = in approximative 99%% confidence limit\n");
81 if (lowcus != -1) { /* single locus curve */
82 fprintf(outfile, "\n\nLocus %li\n(x-axis= Theta\n", lowcus + 1);
83 fprintf(outfile, " y-axis = Recombination-rate,\nunits = log10)\n");
84 for (i = x+1; i >= 0; i--)
85 fprintf(outfile, "%42.42s\n", plane[lowcus][i]);
86 fprintf(outfile, "%42.42s\n", plane[lowcus][x+2]);
87 PAGEFEED;
88 } else {
89 fprintf(outfile, "\nOver all loci\n\n");
90 lowcus = loci;
91 fprintf(outfile, "(x-axis= recombination-rate,");
92 fprintf(outfile, " y-axis = Theta, units = log10)\n");
93 for (i = x; i >= 0; i--)
94 fprintf(outfile, "%42.42s\n", plane[lowcus][i]);
95 PAGEFEED;
96 }
97 } /* plot_surface */
98
99 #if 0
100 void create_loci_plot(world_fmt * world, char **plane, timearchive_fmt * atl, long loci)
101 {
102 long intervals = PLANESIZE;
103 long i, g = 1;
104 nr_fmt *nr;
105 double **pl;
106 double contours[CONTOURLEVELS] = CONTOURS_LOCI;
107
108 nr = (nr_fmt *) calloc(1, sizeof(nr_fmt));
109 for (i = 1; i < loci + 1; i++) {
110 if (g < atl[i].T)
111 g = atl[i].T;
112 }
113 create_nr(nr, atl[loci].numpop, g);
114 pl = (double **) calloc(1, sizeof(double *) * intervals);
115 for (i = 0; i < intervals; i++) {
116 pl[i] = (double *) calloc(1, sizeof(double) * 2 * intervals);
117 }
118 calc_loci_plane(world, nr, atl, pl, loci, contours);
119 fill_plotplane(plane, pl, contours);
120 destroy_nr(nr);
121 print_mathematica(world, pl, intervals, intervals);
122 for (i = 0; i < intervals; i++) {
123 free(pl[i]);
124 }
125 free(pl);
126 }
127 #endif
128
129
create_locus_plot(option_struct * op,data_fmt * data,long lowcus,char ** plane)130 void create_locus_plot(option_struct *op, data_fmt *data, long lowcus, char **plane)
131 {
132 long intervals = PLANESIZE;
133 long i;
134 double **pl;
135 double contours[CONTOURLEVELS] = CONTOURS_LOCUS;
136
137 pl = (double **) calloc(1, sizeof(double *) * intervals);
138 for (i = 0; i < intervals; i++) {
139 pl[i] = (double *) calloc(1, sizeof(double) * 2 * intervals);
140 }
141 calc_locus_plane(op,data,lowcus,pl,contours);
142 fill_plotplane(plane, pl, contours);
143 print_mathematica(pl, intervals, intervals);
144 for (i = 0; i < intervals; i++) {
145 free(pl[i]);
146 }
147 free(pl);
148 }
149
150
151 /*private functions========================================== */
152 /* ---------------------------------------------------------
153 creates memory for archive of timelists for each locus
154 the "locus 0" is reserved for for the current locus, this
155 fake locus is use to speed up the NR-estimation for all chains
156 whereas the loci 1 .. n are used to save the results of the last
157 chain for the combined estimate at the end of the program */
158
create_plotplane(option_struct * op,data_fmt * data,char **** plane)159 void create_plotplane(option_struct *op, data_fmt *data, char**** plane)
160 {
161 short locus, i;
162 long numloci;
163
164 numloci = getdata_numloci(op,data);
165 (*plane) = (char ***) calloc(1, sizeof(char **) * (numloci + 1));
166 for (locus = 0; locus < numloci + 1; locus++) {
167 (*plane)[locus] = (char **) calloc(1, sizeof(char *)
168 * (PLANESIZE + 3));
169 for (i = 0; i < PLANESIZE + 3; i++) {
170 (*plane)[locus][i] = (char *) calloc(1, sizeof(char)
171 * (PLANESIZE + PLANESIZE + 20));
172 }
173 }
174 }
175
176 #if 0
177 void calc_loci_plane(world_fmt * world, nr_fmt * nr, timearchive_fmt * atl, double **pl, long loci, double contours[])
178 {
179 long i, j;
180 double max1 = -DBL_MAX;
181 double max2 = -DBL_MAX;
182 double values[PLANESIZE] = PLANETICKS;
183 long intervals = PLANESIZE;
184 if (world->options->gamma) {
185 nr->param[4] = world->atl[loci + 1].param[4];
186 }
187 for (i = 0; i < intervals; i++) {
188 nr->param[0] = values[i];
189 if (world->options->gamma)
190 calc_gamma(nr);
191 for (j = 0; j < intervals; j++) {
192 nr->param[0] = values[i];
193 nr->param[1] = world->atl[loci + 1].param[1];
194 nr->param[2] = values[j] / values[i];
195 nr->param[3] = world->atl[loci + 1].param[3];
196 calc_loci_like(nr, atl, loci, world->options->gamma);
197 pl[i][j] = nr->llike;
198 if (max1 < nr->llike)
199 max1 = nr->llike;
200 }
201 }
202 nr->param[0] = world->atl[loci + 1].param[0];
203 if (world->options->gamma)
204 calc_gamma(nr);
205 for (i = 0; i < intervals; i++) {
206 for (j = 0; j < intervals; j++) {
207 nr->param[0] = world->atl[loci + 1].param[0];
208 nr->param[1] = values[i];
209 nr->param[2] = world->atl[loci + 1].param[2];
210 nr->param[3] = values[j] / values[i];
211 calc_loci_like(nr, atl, loci, world->options->gamma);
212 pl[i][j + intervals] = nr->llike;
213 if (max2 < nr->llike)
214 max2 = nr->llike;
215 }
216 }
217 contours[0] += max1;
218 contours[1] += max1;
219 contours[2] += max1;
220 contours[3] += max1;
221 contours[4] += max2;
222 contours[5] += max2;
223 contours[6] += max2;
224 contours[7] += max2;
225 }
226 #endif
227
calc_locus_plane(option_struct * op,data_fmt * data,long lowcus,double ** pl,double contours[])228 void calc_locus_plane(option_struct *op, data_fmt *data, long lowcus, double **pl,
229 double contours[])
230 {
231 long intervals = PLANESIZE;
232 long lastchain = totchains - 1;
233 long i, j, chaintype, numloci;
234 double max1 = -DBL_MAX;
235 double max2 = -DBL_MAX;
236 double values[PLANESIZE] = PLANETICKS;
237 double **llike0;
238
239 chaintype = TYPE_CHAIN(lastchain);
240 numloci = getdata_numloci(op,data);
241
242 llike0 = (double **)calloc(1,numloci*sizeof(double *));
243 llike0[0] = (double *)
244 calloc(1,numloci*op->numout[chaintype]*sizeof(double));
245 for(i = 1; i < numloci; i++)
246 llike0[i] = llike0[0] + i*op->numout[chaintype];
247
248 precalc_llike0(op,data,lowcus,lastchain,llike0);
249
250 for (i = 0; i < intervals; i++) {
251 for (j = 0; j < intervals; j++) {
252 pl[i][j] =
253 fmodel_likelihood(op,data,values[j],values[i],llike0,
254 lowcus,lastchain);
255 if (max1 < pl[i][j])
256 max1 = pl[i][j];
257 }
258 }
259 #if 0
260 /* this stuff is for a second population */
261 for (i = 0; i < intervals; i++) {
262 for (j = 0; j < intervals; j++) {
263 nr->param[0] = param[0];
264 nr->param[1] = values[i];
265 nr->param[2] = param[2];
266 nr->param[3] = values[j] / values[i];
267 calc_like(nr, tl, G);
268 pl[i][j + intervals] = nr->llike;
269 if (max2 < nr->llike)
270 max2 = nr->llike;
271 }
272 }
273 #endif
274 contours[0] += max1;
275 contours[1] += max1;
276 contours[2] += max1;
277 contours[3] += max1;
278 contours[4] += max2;
279 contours[5] += max2;
280 contours[6] += max2;
281 contours[7] += max2;
282
283 free(llike0[0]);
284 free(llike0);
285
286 } /* calc_locus_plane */
287
fill_plotplane(char ** plane,double ** pl,double contours[])288 void fill_plotplane(char **plane, double **pl, double contours[])
289 {
290 long i, j, zz = 0;
291
292 char line[100];
293 long val[PLANEBIGTICKS] = PLANEBIGTICKVALUES;
294 long intervals = PLANESIZE;
295 for (i = 0; i < intervals; i++) {
296 if (i % 7)
297 line[i] = '-';
298 else
299 line[i] = '+';
300 }
301 line[i] = '\0';
302 sprintf(plane[0], " +%s+ +%s+ ", line, line);
303 for (i = 0; i < intervals; i++) {
304 memset(plane[i + 1], ' ', sizeof(char) * (intervals + intervals + 20));
305 plane[i + 1][intervals + intervals + 19] = '\0';
306 if (!((i) % 7)) {
307 sprintf(plane[i + 1], " %2.0li +", val[zz++]);
308 if(val[zz-1]==0)
309 plane[i+1][3]='0';
310 }
311 else
312 plane[i + 1][5] = '|';
313 for (j = 0; j < intervals; j++) {
314 if (pl[i][j] < contours[1]) {
315 if (pl[i][j] < contours[2]) {
316 if (pl[i][j] < contours[3]) {
317 plane[i + 1][j + 6] = ' ';
318 }
319 else {
320 plane[i + 1][j + 6] = '-';
321 }
322 }
323 else {
324 plane[i + 1][j + 6] = '+';
325 }
326 }
327 else {
328 if (pl[i][j] < contours[0] - EPSILON)
329 plane[i + 1][j + 6] = '*';
330 else {
331 plane[i + 1][j + 6] = 'X';
332 }
333 }
334 }
335 if ((i) % 7) {
336 plane[i + 1][j + 6] = '|';
337 plane[i + 1][j + 11] = '|';
338 }
339 else {
340 plane[i + 1][j + 6] = '+';
341 plane[i + 1][j + 11] = '+';
342 }
343 #if 0 /* this is stuff for the second population */
344 for (j = intervals + 7 + 5 /*- 1*/; j < intervals + 7 + 5 + intervals /*- 2*/; j++) {
345 z = j - 7 - 5 /*+ 1*/;
346 if (pl[i][z] < contours[5]) {
347 if (pl[i][z] < contours[6]) {
348 if (pl[i][z] < contours[7]) {
349 plane[i + 1][j] = ' ';
350 }
351 else {
352 plane[i + 1][j] = '-';
353 }
354 }
355 else {
356 plane[i + 1][j] = '+';
357 }
358 }
359 else {
360 if (pl[i][z] < contours[4] - EPSILON)
361 plane[i + 1][j] = '*';
362 else {
363 plane[i + 1][j] = 'X';
364 }
365 }
366 }
367 if ((i) % 7) {
368 plane[i + 1][j] = '|';
369 plane[i + 1][j + 1] = '\0';
370 }
371 else {
372 plane[i + 1][j] = '+';
373 plane[i + 1][j + 1] = '\0';
374 }
375 #endif
376 }
377 sprintf(plane[intervals+1], " +%s+ +%s+", line, line);
378 sprintf(plane[intervals+2], " -3 -2 -1 0 1 2 -3 -2 -1 0 1 2");
379 }
380
print_mathematica(double ** plane,long x,long y)381 void print_mathematica(double **plane, long x, long y)
382 {
383 long i, j;
384 static long number = 1;
385 FILE *mathfile;
386
387 mathfile = fopen("mathfile","w+");
388
389 fprintf(mathfile, "\n\nlocus%-li={{\n", number++);
390 for (i = 0; i < x-1; i++) {
391 fprintf(mathfile, "{");
392 for (j = 0; j < y-1; j++) {
393 fprintf(mathfile, "%20.20g,", plane[i][j]);
394 }
395 fprintf(mathfile, "%20.20g},\n", plane[i][j]);
396 }
397 fprintf(mathfile, "{");
398 for (j = 0; j < y-1; j++) {
399 fprintf(mathfile, "%20.20g,", plane[i][j]);
400 }
401 fprintf(mathfile, "%20.20g}},{\n", plane[i][j]);
402 for (i = 0; i < x-1; i++) {
403 fprintf(mathfile, "{");
404 for (j = y; j < 2 * y-1; j++) {
405 fprintf(mathfile, "%20.20g,", plane[i][j]);
406 }
407 fprintf(mathfile, "%20.20g},\n", plane[i][j]);
408 }
409 fprintf(mathfile, "{");
410 for (j = y; j < 2 * y-1; j++) {
411 fprintf(mathfile, "%20.20g,", plane[i][j]);
412 }
413 fprintf(mathfile, "%20.20g}}}\n", plane[i][j]);
414
415 fclose(mathfile);
416
417 } /* print_mathematica */
418
419
free_plotplane(option_struct * op,data_fmt * data,char *** plane)420 void free_plotplane(option_struct *op, data_fmt *data, char ***plane)
421 {
422 long locus, i, numloci;
423
424 numloci = getdata_numloci(op,data);
425
426 for (locus = 0; locus < numloci + 1; locus++) {
427 for (i = 0; i < PLANESIZE + 3; i++) free(plane[locus][i]);
428 free(plane[locus]);
429 }
430
431 free(plane);
432
433 } /* free_plotplane */
434
435
print_locusplot(option_struct * op,data_fmt * data,long lowcus)436 void print_locusplot (option_struct *op, data_fmt *data, long lowcus)
437 {
438 char ***plane;
439 long numloci;
440
441 numloci = getdata_numloci(op,data);
442
443 create_plotplane(op,data,&plane);
444 if (lowcus!= -1)
445 create_locus_plot(op,data,lowcus,plane[lowcus]);
446 else create_locus_plot(op,data,lowcus,plane[numloci]);
447 plot_surface(op,data,lowcus,plane,(long)PLANESIZE);
448 free_plotplane(op,data,plane);
449
450 } /* print_locusplot */
451