1 /*------------------------------------------------------
2 Maximum likelihood estimation
3 of migration rate and effectice population size
4 using a Metropolis-Hastings Monte Carlo algorithm
5 -------------------------------------------------------
6 mighistogram R O U T I N E S that use less memory
7
8 Peter Beerli 2006, Tallahassee
9 beerli@fsu.edu
10
11 Copyright 2006 Peter Beerli, Tallahassee
12
13 Permission is hereby granted, free of charge, to any person obtaining
14 a copy of this software and associated documentation files (the
15 "Software"), to deal in the Software without restriction, including
16 without limitation the rights to use, copy, modify, merge, publish,
17 distribute, sublicense, and/or sell copies of the Software, and to
18 permit persons to whom the Software is furnished to do so, subject
19 to the following conditions:
20
21 The above copyright notice and this permission notice shall be
22 included in all copies or substantial portions of the Software.
23
24 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
25 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
26 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
27 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
28 ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
29 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
30 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
31
32
33 $Id$
34 */
35 /*! \file migevents.c
36
37 this file contains functions that repport the numbers and times of coalescence and migration events
38 Results will be printed as a histogram over time and also will print the histogram to file and
39 print a summary table
40 */
41 #include <stdlib.h>
42
43 #include "migevents.h"
44 #include "skyline.h"
45 #include "random.h"
46 #include "tools.h"
47 #include "sighandler.h"
48 #include "world.h"
49 #ifdef PRETTY
50 #include "pretty.h"
51 #endif
52 #ifdef MPI
53 #include "migrate_mpi.h"
54 #else
55 extern int myID;
56 #endif
57
58 void setup_mig_coal_events (world_fmt * world, option_fmt * options);
59
60 ///
61 /// allocates the necessary memory for recording of the events through time
62 void
setup_mighist(world_fmt * world,option_fmt * options)63 setup_mighist (world_fmt * world, option_fmt * options)
64 {
65 if (world->options->mighist)
66 {
67 world->mighistloci = (mighistloci_fmt *)
68 mycalloc (world->loci+1, sizeof (mighistloci_fmt));
69 world->mighistlocinum = 0;
70 setup_mig_coal_events(world,options);
71 }
72 else
73 world->mighistloci = NULL;
74 }
75
76 ///
77 /// destructor of the container that records events through time
78 void
destroy_mighist(world_fmt * world)79 destroy_mighist (world_fmt * world)
80 {
81 long locus;
82 long pop;
83 if (world->options->mighist)
84 {
85 for(locus=0; locus < world->loci; locus++)
86 {
87 for(pop=0;pop < world->numpop2; pop++)
88 {
89 myfree(world->mighistloci[locus].migeventbins[pop]);
90 }
91 myfree(world->mighistloci[locus].migeventbins);
92 // if(world->options->skyline)
93 // myfree(world->mighistloci[locus].eventbins);
94 }
95 myfree(world->mighistloci);
96 }
97 }
98
99 void
minmax(histogram_fmt * hist,float * tempmin,float * tempmax)100 minmax (histogram_fmt * hist, float *tempmin, float *tempmax)
101 {
102 long i;
103 MYREAL tmp1, tmp2;
104 MYREAL tmpmin = MYREAL_MAX;
105 MYREAL tmpmax = -MYREAL_MAX;
106
107 for (i = 0; i < hist->count; i++)
108 {
109 if ((tmp1 = hist->time[i]) < tmpmin)
110 tmpmin = tmp1;
111 if ((tmp2 = hist->time[i]) > tmpmax)
112 tmpmax = tmp2;
113 }
114 *tempmax = (float) tmpmax;
115 *tempmin = (float) tmpmin;
116 }
117
118 ///
119 /// print and calculate means, medians, of migration and coalescence event times
120 /// over all times and of the last (the MRCA) event.
121 void
print_mighist_output(FILE * out,world_fmt * world,MYREAL * sums,boolean mrca)122 print_mighist_output (FILE * out, world_fmt * world, MYREAL *sums, boolean mrca)
123 {
124 const long column = (mrca ? 1 : 0 );
125 const long numpop2 = world->numpop2;
126 const long startpop = (world->options->mighist_all ? 0 : world->numpop);
127 const long endpop = (mrca ? world->numpop : numpop2);
128 const long loci1 = ((world->loci == 1) ? 1 : world->loci + 1);
129 const long loci1pop = loci1 * numpop2;
130 const long sumloc = ((world->loci == 1) ? 0 : world->loci);
131
132 long locus;
133 long pop;
134 long i;
135 long lp;
136 long frompop;
137 long topop;
138 long * eventbinnum = NULL;
139 long * eventbinmax = NULL;
140
141 duo ** eventbins;
142
143 float *total;
144 float *meantime;
145 float *mediantime;
146 float *stdtime;
147 float *freq;
148 float *values;
149 float **allvalues=NULL;
150 float ntotal;
151 float eventbinsize;
152 float age;
153 float binfreq;
154
155 eventbinmax = (long *) mycalloc(numpop2, sizeof(long));
156
157
158 if(loci1>1)
159 {
160 for (pop = startpop; pop < endpop; pop++)
161 {
162 for (locus = 0; locus < world->loci; locus++)
163 {
164 eventbinnum = world->mighistloci[locus].migeventbinnum;
165 if(eventbinnum[pop] > eventbinmax[pop])
166 eventbinmax[pop] = eventbinnum[pop];
167 }
168 }
169 allvalues = (float **) mycalloc(world->numpop2,sizeof(float*));
170 for (pop = startpop; pop < endpop; pop++)
171 {
172 allvalues[pop] = (float *) mycalloc(eventbinmax[pop]+1,sizeof(float));
173 }
174 }
175 // initialize the vectors used for the table
176 values = (float *) mycalloc(1, sizeof(float));
177 total = (float *) mycalloc(loci1, sizeof(float));
178 meantime = (float *) mycalloc(loci1pop * 4, sizeof(float));
179 mediantime = meantime + (loci1pop);
180 stdtime = mediantime + (loci1pop);
181 freq = stdtime + (loci1pop);
182
183 for (locus = 0; locus < world->loci; locus++)
184 {
185 eventbins = world->mighistloci[locus].migeventbins;
186 eventbinnum = world->mighistloci[locus].migeventbinnum;
187 eventbinsize = world->mighistloci[locus].eventbinsize;
188
189 for (pop = startpop; pop < endpop; pop++)
190 {
191 age = eventbinsize / 2.;
192 values = (float *) myrealloc(values, eventbinnum[pop] * sizeof(float));
193
194 lp = locus * numpop2 + pop;
195
196 meantime[lp] = 0.f;
197 stdtime[lp] = 0.f;
198 freq[lp] = 0.f;
199
200 for (i = 0; i < eventbinnum[pop]; i++)
201 {
202 binfreq = (float) (eventbins[pop][i][column]);
203 freq[lp] += binfreq; // total frequency of events per parameter
204 total[locus] += binfreq; // total of frequency this should add up to world->loci
205 values[i] = binfreq * age; // frequency times age [how many times was this age recorded.]
206 meantime[lp] += values[i];
207 stdtime[lp] += age * age * binfreq;
208 age += eventbinsize;
209 if(loci1 > 1)
210 {
211 allvalues[pop][i] += values[i];
212 }
213 if(freq[lp] < 0.499999)
214 {
215 mediantime[lp] = age;
216 }
217 }
218 // freq[lp] should be one, but the divisions are done to fix accumulated rounding issues.
219 if(freq[lp]>0.0)
220 {
221 meantime[lp] /= freq[lp];
222 stdtime[lp] = ((stdtime[lp]/freq[lp]) - (meantime[lp]*meantime[lp]));
223 }
224 if(loci1>1)
225 {
226 meantime[sumloc * world->numpop2 + pop] += meantime[lp]/world->loci;
227 stdtime[sumloc * world->numpop2 + pop] += stdtime[lp]/world->loci;
228 mediantime[sumloc * world->numpop2 + pop] += mediantime[lp]/world->loci;
229 }
230 }
231 // calculating the frequences over all parameters using the sums
232 ntotal = 0.;
233 for (pop = startpop; pop < endpop; pop++)
234 {
235 lp = locus * numpop2 + pop;
236 ntotal += sums[lp];
237 freq[lp] = sums[lp];
238 stdtime[lp] = sqrt(stdtime[lp]);
239 }
240 for (pop = startpop; pop < endpop; pop++)
241 {
242 lp = locus * numpop2 + pop;
243 freq[lp] /= ntotal;
244 if(loci1>1)
245 freq[sumloc * numpop2 + pop] += freq[lp]/world->loci;
246 }
247 if(loci1>1)
248 total[sumloc] += total[locus];
249 }
250
251 myfree(values);
252
253 if(loci1 > 1)
254 {
255 for (pop = startpop; pop < endpop; pop++)
256 {
257 lp = sumloc * numpop2 + pop;
258 //mediantime[lp] = -1.; //DEBUG: allvalues[pop][(long) (eventbinnum[pop]/2)]/freq[lp];
259 stdtime[lp] = sqrt(stdtime[lp]);
260 }
261 }
262
263 if(mrca)
264 {
265 FPRINTF (out, "\n\nTime and probability of most recent common ancestor\n");
266 FPRINTF (out, "===================================================\n\n");
267 }
268 else
269 {
270 FPRINTF (out, "\n\nApproximate Summary of %s events\n",
271 world->options->mighist_all ? "coalescence and migration" : "migration");
272 FPRINTF (out, "===============%s================\n\n",
273 world->options->mighist_all ? "=========================" : "=========");
274 }
275 for (locus = 0; locus < loci1; locus++)
276 { /* Each locus + Summary */
277 if (locus != world->loci)
278 FPRINTF (out, "Locus %li\n", locus + 1);
279 else
280 FPRINTF (out, "Over all loci\n");
281 FPRINTF (out,
282 "---------------------------------------------------------\n");
283 FPRINTF (out,
284 "Population Time Frequency\n");
285 FPRINTF (out, " -----------------------------\n");
286 FPRINTF (out, "From To Average Median Std\n");
287 FPRINTF (out,
288 "---------------------------------------------------------\n");
289 for (pop = (world->options->mighist_all ? 0 : world->numpop); pop < (mrca ? world->numpop : world->numpop2); pop++)
290 {
291 if(world->bayes->map[pop][1] == INVALID)
292 continue;
293 lp = locus * world->numpop2 + pop;
294 m2mm(pop,world->numpop,&frompop,&topop);
295 if(freq[lp]<=0.0)
296 {
297 FPRINTF (out,
298 "%4li %4li ------- ------- ------- %3.5f\n",
299 frompop + 1,
300 topop + 1,
301 0.0);
302 }
303 else
304 {
305 FPRINTF (out,
306 "%4li %4li %3.5f %3.5f %3.5f %3.5f\n",
307 frompop + 1,
308 topop + 1,
309 meantime[lp], mediantime[lp],
310 stdtime[lp],
311 freq[lp]);
312 }
313 }
314 FPRINTF (out,
315 "---------------------------------------------------------\n");
316 FPRINTF (out, "\n");
317 }
318 #ifdef PRETTY
319 pdf_print_time_table (world, meantime, mediantime, stdtime, freq, mrca);
320 #endif
321
322 myfree(total);
323 myfree(meantime);
324 myfree(allvalues);
325 myfree(eventbinmax);
326 }
327
328
329
330 ///
331 /// gather the events from the genealogy for each timeinterval.
332 /// This function constructs a count-histogram for all event types
333 /// and a specific histogram for the last coalescent (the MRCA)
calculate_event_values(duo ** eventbins,long * eventbinnum,MYREAL eventinterval,MYREAL interval,MYREAL age,long from,long to,long * lineages,long numpop,boolean is_last)334 void calculate_event_values(duo **eventbins, long *eventbinnum, MYREAL eventinterval,
335 MYREAL interval, MYREAL age, long from, long to,
336 long * lineages, long numpop, boolean is_last)
337 {
338 MYREAL inv_eventinterval = 1. / eventinterval;
339 long bin = (long) (age * inv_eventinterval); // this bin
340 // long k;
341 long i;
342 long pop;
343 MYREAL val;
344 long lastbin = (long) ((age - interval) * inv_eventinterval);
345 long diff;
346 duo *bins=NULL;
347 long allocsize=1;
348 if(lastbin < 0)
349 {
350 lastbin=0;
351 }
352 // if(bin > 10000 || bin < 0) // tooo many bins and it overflows probably
353 // //I should stop with some bins like 10000
354 // {
355 // warning("Bin (%li) is smaller than zero or bigger than 10000", bin);
356 // return;
357 // }
358 diff = bin - lastbin;
359 val = 1.0;
360 if(from == to)
361 {
362 pop = to;
363 //fprintf(stdout,"%i> c bin=%li ebin[%li]=%li\n",myID, bin, pop, eventbinnum[pop]);
364 }
365 else
366 {
367 pop = mm2m(from, to, numpop);
368 //fprintf(stdout,"%i> m bin=%li ebin[%li]=%li\n",myID, bin, pop, eventbinnum[pop]);
369 }
370 //fflush(stdout);
371 if(bin>=eventbinnum[pop])
372 {
373 allocsize = bin+1;
374 eventbins[pop] = (duo *) myrealloc(eventbins[pop], allocsize * sizeof(duo));
375 // fprintf(stdout,"%i> ral: as=%li, old=%li, last val=(%f,%f)\n",myID,allocsize, eventbinnum[pop],
376 // eventbins[pop][eventbinnum[pop]-1][0],eventbins[pop][eventbinnum[pop]-1][1]);
377 for(i=eventbinnum[pop];i < allocsize; i++)
378 {
379 eventbins[pop][i][0] = 0.;
380 eventbins[pop][i][1] = 0.;
381 }
382 eventbinnum[pop] = allocsize;
383 }
384 bins = eventbins[pop];
385 // fprintf(stdout,"m diff=%li ebin[pop]=%li, from=%i to=%li val=%f interval=%f bin=%li\n",diff , eventbinnum[pop], from, pop, val, interval, bin);
386
387 if(diff < 0)
388 error("time difference is negative -- not possible");
389 bins[bin][0] += val;
390 if(is_last)
391 bins[bin][1] += val;
392 }
393
394 ///
395 /// set up the event plot histogram containers, only when also the migration histograms are recorded
396 /// this function needs to be called by setup_mighist() function
397 void
setup_mig_coal_events(world_fmt * world,option_fmt * options)398 setup_mig_coal_events (world_fmt * world, option_fmt * options)
399 {
400 long locus, i, j;
401 long allocsize = 1;
402 MYREAL binsize = world->options->eventbinsize; // in mutational units (=time scale)
403 if (world->options->mighist)
404 {
405 for (locus = 0; locus < world->loci; locus++)
406 {
407 world->mighistloci[locus].eventbinsize = binsize;
408 if(world->mighistloci[locus].migeventbins == NULL)
409 world->mighistloci[locus].migeventbins = (duo **) mycalloc (world->numpop2, sizeof (duo *));
410 if(world->mighistloci[locus].migeventbinnum == NULL)
411 world->mighistloci[locus].migeventbinnum = (long *) mycalloc (world->numpop2, sizeof (long));
412
413 for (i = 0; i < world->numpop2; i++)
414 {
415 world->mighistloci[locus].migeventbinnum[i] = allocsize;
416 if(world->mighistloci[locus].migeventbins[i] == NULL)
417 world->mighistloci[locus].migeventbins[i] = (duo *) mycalloc (allocsize, sizeof (duo));
418 for(j=0;j<allocsize;j++)
419 {
420 world->mighistloci[locus].migeventbins[i][j][0] = 0.;
421 world->mighistloci[locus].migeventbins[i][j][1] = 0.;
422 }
423 }
424 }
425 }
426 }
427
428 ///
429 /// Destroy the event plot histogram container
430 void
destroy_mig_coal_events(world_fmt * world)431 destroy_mig_coal_events (world_fmt * world)
432 {
433 long locus, i;
434 if (world->options->mighist)
435 {
436 for (locus = 0; locus < world->loci; locus++)
437 {
438 for (i = 0; i < world->numpop2; i++)
439 {
440 myfree(world->mighistloci[locus].migeventbins[i]);
441 }
442 myfree(world->mighistloci[locus].migeventbins);
443 myfree(world->mighistloci[locus].migeventbinnum);
444 }
445 }
446 }
447
448
print_event_values_list(FILE * file,long locus,duo ** eventbins,MYREAL eventbinsize,long * eventbinnum,long numpop)449 void print_event_values_list(FILE *file, long locus, duo **eventbins, MYREAL eventbinsize, long *eventbinnum, long numpop)
450 {
451 long i;
452 long pop;
453 long frompop;
454 long topop;
455 long numpop2 = numpop * numpop;
456 MYREAL age;
457
458 for(pop = 0; pop < numpop2; pop++)
459 {
460 //xcode age = 0.;
461 if(pop < numpop)
462 {
463 fprintf(file,"\nLocus: %li Parameter: %s_%li\n", locus+1, "Theta",pop+1);
464 }
465 else
466 {
467 m2mm(pop,numpop,&frompop,&topop);
468 fprintf(file,"\nLocus: %li Parameter: %s_(%li,%li)\n", locus+1, "M", frompop+1, topop+1);
469 }
470 fprintf(file,"Time Frequency of visit Frequency of MRCA\n");
471 fprintf(file,"--------------------------------------------------\n");
472 age = eventbinsize/2.;
473 for(i = 0; i < eventbinnum[pop]; i++)
474 {
475 fprintf(file,"%10.10f %10.10f %10.10f\n", age, eventbins[pop][i][0], eventbins[pop][i][1]);
476 age += eventbinsize;
477 }
478 }
479 }
480
481 ///
482 /// read events from mighistfile and push them into the mighistlocus structure
483 #define SOME_ELEMENTS 100
read_event_values_fromfile(FILE * file,world_fmt * world)484 void read_event_values_fromfile(FILE *file, world_fmt *world)
485 {
486 long locus;
487 long pop;
488 long i;
489 // MYREAL age;
490 char *input;
491 char *inptr;
492 duo **eventbins;
493 long *eventbinnum;
494 //MYREAL eventbinsize;
495 long spacer=0;
496 long *allocsize;
497 input = (char *) mycalloc(SUPERLINESIZE, sizeof(char));
498 allocsize = (long *) mycalloc(world->numpop2*world->loci,sizeof(long));
499 for(i=0;i<world->numpop2*world->loci;i++)
500 {
501 allocsize[i]=1;
502 }
503 while(FGETS(input,SUPERLINESIZE,file) != EOF)
504 {
505 // grab the commentlines
506 while(input[0]=='#')
507 {
508 FGETS(input,LINESIZE,file);
509 }
510 // read the mighistfile
511 if(input !=NULL)
512 {
513 inptr = input;
514 locus = atol(strsep(&inptr,"\t"))-1;
515 spacer = locus * world->numpop2;
516 if(locus == -1)
517 error("help");
518 pop = atol(strsep(&inptr,"\t"))-1;
519 i = atol(strsep(&inptr,"\t"))-1;
520 //age = atof(strsep(&inptr,"\t"));
521 (void) strsep(&inptr,"\t");
522 eventbins = world->mighistloci[locus].migeventbins;
523 eventbinnum = world->mighistloci[locus].migeventbinnum;
524 //xcode eventbinsize = world->mighistloci[locus].eventbinsize;
525
526 //xcode if(i==0)
527 //xcode {
528 //xcode eventbinsize = 2 * age;
529 //xcode }
530 eventbinnum[pop] = i+1;
531 if(allocsize[spacer + pop] < eventbinnum[pop])
532 {
533 allocsize[spacer + pop] += SOME_ELEMENTS;
534 eventbins[pop] = (duo *) myrealloc(eventbins[pop],allocsize[spacer + pop] * sizeof(duo));
535 }
536 eventbins[pop][i][0] = atof(strsep(&inptr,"\t"));
537 eventbins[pop][i][1] = atof(strsep(&inptr,"\t"));
538 }
539 }
540 myfree(allocsize);
541 }
542
543
print_event_values_tofile(FILE * file,world_fmt * world)544 void print_event_values_tofile(FILE *file, world_fmt *world)
545 {
546 long i;
547 long pop;
548 long frompop;
549 long topop;
550 long locus;
551 long sumloc;
552 long numpop = world->numpop;
553 long numpop2 = world->numpop2;
554 MYREAL age;
555 duo **eventbins;
556 long *eventbinnum;
557 MYREAL eventbinsize;
558
559 fprintf(file,"# Raw record of the events histogram for all parameters and all loci\n");
560 fprintf(file,"# The time interval is set to %f\n", world->options->eventbinsize);
561 fprintf(file,"# produced by the program %s (http://popgen.scs.fsu.edu/migrate.hml)\n",
562 MIGRATEVERSION);
563 fprintf(file,"# written by Peter Beerli 2006-2007, Tallahassee,\n");
564 fprintf(file,"# if you have problems with this file please email to beerli@fsu.edu\n");
565 fprintf(file,"#\n");
566 fprintf(file,"# Order of the parameters:\n");
567 fprintf(file,"# Parameter-number Parameter\n");
568 for(pop=0;pop<numpop2;pop++)
569 {
570 if(pop < numpop)
571 {
572 fprintf(file,"# %6li %s_%li\n", pop+1, "Theta",pop+1);
573 }
574 else
575 {
576 m2mm(pop,numpop,&frompop,&topop);
577 fprintf(file,"# %6li %s_(%li,%li)\n", pop+1, (world->options->usem ? "M" : "xNm"), frompop+1, topop+1);
578 }
579 }
580 fprintf(file,"#\n#----------------------------------------------------------------------------\n");
581 fprintf(file,"# Locus Parameter-number Bin Age Events-frequency(*) MRCA-frequency(*)\n");
582 fprintf(file,"#----------------------------------------------------------------------------\n");
583 fprintf(file,"# (*) values with -1 were NEVER visited\n");
584 fprintf(file,"# @@@@@Locus param bin age freq mrcafreq\n");
585 if(world->loci>1)
586 {
587 sumloc =1;
588 fprintf(file,"# Locus %li is sum over all loci, when there are more than 1 locus\n", world->loci+1);
589 }
590 else
591 {
592 sumloc = 0;
593 }
594
595 for(locus=0; locus < world->loci + sumloc; locus++)
596 {
597 if(!world->data->skiploci[locus])
598 {
599 eventbins = world->mighistloci[locus].migeventbins;
600 eventbinnum = world->mighistloci[locus].migeventbinnum;
601 eventbinsize = world->mighistloci[locus].eventbinsize;
602 for(pop = 0; pop < numpop2; pop++)
603 {
604 age = eventbinsize / 2.;
605 for(i = 0; i < eventbinnum[pop]; i++)
606 {
607 fprintf(file,"%li\t%li\t%li\t%10.10f\t%10.10f\t%10.10f\n",
608 locus+1, pop+1, i+1, age, eventbins[pop][i][0], eventbins[pop][i][1]);
609 age += eventbinsize;
610 }
611 }
612 }
613 }
614 }
615
616 ///
617 /// calculates frequencies for all event bins (over all and MRCA)
prepare_event_values(world_fmt * world,MYREAL * sums0,MYREAL * sums1)618 void prepare_event_values(world_fmt *world, MYREAL *sums0,MYREAL *sums1)
619 {
620 MYREAL sum0;
621 MYREAL sum1;
622 long z = 0;
623 long locus;
624 long pop;
625 long i;
626 long numpop2 = world->numpop2;
627 duo **eventbins;
628 duo **eventbins_all;
629 long *eventbinnum;
630 long eventbinnum_allmax=0;
631
632 for(locus=0; locus < world->loci; locus++)
633 {
634 if(!world->data->skiploci[locus])
635 {
636 eventbins = world->mighistloci[locus].migeventbins;
637 eventbinnum = world->mighistloci[locus].migeventbinnum;
638 for(pop = 0; pop < numpop2; pop++)
639 {
640 sum0 = (MYREAL) 0.0;
641 sum1 = (MYREAL) 0.0;
642
643 if(eventbinnum[pop] > eventbinnum_allmax)
644 eventbinnum_allmax = eventbinnum[pop];
645
646 for(i = 0; i < eventbinnum[pop]; i++)
647 {
648 if(eventbins[pop][i][0] <= 0.0)
649 {
650 continue;
651 }
652 else
653 {
654 sum0 += eventbins[pop][i][0];
655 }
656
657 if(eventbins[pop][i][1] <= 0.0)
658 {
659 continue;
660 }
661 else
662 {
663 sum1 += eventbins[pop][i][1];
664 }
665 }
666 // calculate frequency from bincount data
667 for(i = 0; i < eventbinnum[pop]; i++)
668 {
669 if(eventbins[pop][i][0] <= 0.0)
670 continue;
671 eventbins[pop][i][0] /= sum0;
672 if(eventbins[pop][i][1] <= 0.0)
673 continue;
674 eventbins[pop][i][1] /= sum1;
675 }
676 sums0[z]=sum0;
677 sums1[z++]=sum1;
678 }
679 }
680 }
681 // for all loci, even when there is only one
682 if(world->mighistloci[world->loci].migeventbins == NULL)
683 world->mighistloci[world->loci].migeventbins = (duo **) mycalloc(world->numpop2,sizeof(duo *));
684 if(world->mighistloci[world->loci].migeventbinnum == NULL)
685 world->mighistloci[world->loci].migeventbinnum = (long *) mycalloc(world->numpop2,sizeof(long));
686 world->mighistloci[world->loci].eventbinsize = world->mighistloci[0].eventbinsize;
687 for(pop=0; pop< world->numpop2 ; pop++)
688 {
689 if(world->mighistloci[world->loci].migeventbins[pop] == NULL)
690 world->mighistloci[world->loci].migeventbins[pop] = (duo *) mycalloc(eventbinnum_allmax,sizeof(duo));
691 world->mighistloci[world->loci].migeventbinnum[pop] = eventbinnum_allmax;
692 }
693 eventbins_all = world->mighistloci[world->loci].migeventbins;
694 for (locus = 0; locus < world->loci; locus++)
695 {
696 if(!world->data->skiploci[locus])
697 {
698 eventbinnum = world->mighistloci[locus].migeventbinnum;
699 eventbins = world->mighistloci[locus].migeventbins;
700 for(pop=0; pop< world->numpop2 ; pop++)
701 {
702 for(i=0 ; i < eventbinnum[pop]; i++)
703 {
704 if(eventbins[pop][i][0] > 0.0)
705 {
706 eventbins_all[pop][i][0] += eventbins[pop][i][0]/world->loci;
707 }
708 if(eventbins[pop][i][1] > 0.0)
709 {
710 eventbins_all[pop][i][1] += eventbins[pop][i][1]/world->loci;
711 }
712 }
713 }
714 }
715 }
716 }
717
718 ///
719 /// print titles for events over time
print_event_values_title(FILE * file,boolean progress)720 void print_event_values_title(FILE *file, boolean progress)
721 {
722 if(progress)
723 {
724 fprintf(file,"\n\nEvents over time\n");
725 fprintf(file,"--------------------\n\n");
726 }
727 }
728
729 ///
730 /// print event statistics to oufile, mighistfile and also PDF
print_event_values(world_fmt * world)731 void print_event_values(world_fmt * world)
732 {
733 long locus;
734 long sumloc = (world->loci > 1) ? 1 : 0;
735 MYREAL *sums0;
736 MYREAL *sums1;
737
738 if(world->options->mighist)
739 {
740
741 sums0 = (MYREAL *) mycalloc(world->loci * world->numpop2 * 2,sizeof(MYREAL));
742 sums1 = sums0 + (world->loci * world->numpop2);
743
744 if(world->options->datatype != 'g')
745 {
746 // prepare event histogram for printing and get total of observations
747 // for all (sums0) and for the mrcas only (sums1)
748 prepare_event_values(world,sums0, sums1);
749
750 // print event to file
751 print_event_values_tofile(world->mighistfile, world);
752 }
753 else
754 {
755 read_event_values_fromfile(world->mighistfile,world);
756 prepare_event_values(world,sums0, sums1);
757 }
758 // print title to screen
759 print_event_values_title(stdout, world->options->progress);
760
761 // print title to ascii-outfile
762 print_event_values_title(world->outfile, TRUE);
763 #ifdef PRETTY
764 if(world->options->verbose)
765 pdf_print_eventtime_table(world);
766 #endif
767 // print content to stdout and to ascii - outfile [it seems not useful in the output file
768 if(world->options->verbose)
769 {
770 for(locus=0; locus < world->loci+sumloc; locus++)
771 {
772 if(!world->data->skiploci[locus])
773 {
774 print_event_values_list(stdout, locus, world->mighistloci[locus].migeventbins,
775 world->mighistloci[locus].eventbinsize,
776 world->mighistloci[locus].migeventbinnum, world->numpop);
777 print_event_values_list(world->outfile, locus, world->mighistloci[locus].migeventbins,
778 world->mighistloci[locus].eventbinsize,
779 world->mighistloci[locus].migeventbinnum, world->numpop);
780 }
781 }
782 }
783 #ifdef PRETTY
784 pdf_event_histogram(world->loci, world->numpop2, world);
785 #endif
786 // this include also printing to the pretty interface
787 print_mighist_output (world->outfile, world, sums0, FALSE);
788 print_mighist_output (world->outfile, world, sums1, TRUE);
789 myfree(sums0);
790 }
791 }
792
793 /// store the event statistics
794 void
store_events(world_fmt * world,timelist_fmt * ltl,long np,long rep)795 store_events (world_fmt * world, timelist_fmt * ltl, long np, long rep)
796 {
797 long j;
798 long T;
799 MYREAL t;
800 mighistloci_fmt *bb = NULL;
801 if (world->in_last_chain && world->options->mighist)
802 {
803 world->options->mighist_counter = 0;
804
805 if(!world->data->skiploci[world->locus])
806 {
807 bb = &(world->mighistloci[world->locus]);
808
809 if((*ltl).tl[0].eventnode->type != 't')
810 {
811 calculate_event_values(bb->migeventbins, bb->migeventbinnum, bb->eventbinsize,
812 (*ltl).tl[0].age,
813 (*ltl).tl[0].age, (*ltl).tl[0].from, (*ltl).tl[0].to,
814 (*ltl).tl[0].lineages, world->numpop, FALSE);
815 // calculate skyline values
816 if(world->options->skyline)
817 {
818 if((*ltl).tl[0].eventnode->visited)
819 {
820 calculate_expected_values(bb->eventbins, bb->eventbinnum, bb->eventbinsize,
821 (*ltl).tl[0].age,
822 (*ltl).tl[0].age, (*ltl).tl[0].from, (*ltl).tl[0].to,
823 (*ltl).tl[0].lineages, world->numpop, world);
824 }
825 }
826 }
827 T = (*ltl).T - 1;
828 for (j = 1; j < T; j++)
829 {
830 t = (*ltl).tl[j].age - (*ltl).tl[j - 1].age;
831 // if((*ltl).tl[j].eventnode->type=='t')
832 // continue;
833 if((*ltl).tl[j].eventnode->type != 't')
834 {
835 // calculate event values
836 calculate_event_values(bb->migeventbins, bb->migeventbinnum, bb->eventbinsize,
837 t,
838 (*ltl).tl[j].age, (*ltl).tl[j].from, (*ltl).tl[j].to,
839 (*ltl).tl[j].lineages, world->numpop, is_same(j,T-1));
840
841 // calculate skyline values
842 if(world->options->skyline)
843 {
844 if((*ltl).tl[j].eventnode->visited)
845 {
846 calculate_expected_values(bb->eventbins, bb->eventbinnum, bb->eventbinsize, t,
847 (*ltl).tl[j].age, (*ltl).tl[j].from, (*ltl).tl[j].to,
848 (*ltl).tl[j].lineages, world->numpop, world);
849 }
850 }
851 }
852 }
853 }
854 }
855 }
856