1 /*! \file migration.h*/ 2 #ifndef MIGRATION_HEADER /*migration.h */ 3 #define MIGRATION_HEADER 4 5 #include "definitions.h" 6 /*----------------------------------------------------------------- 7 Maximum likelihood estimation of migration rates 8 using coalescent trees 9 10 Peter Beerli 11 School of Computational Science 12 Department of Biological Sciences 13 Florida State University 14 Tallahassee FL 32306-4120 15 beerli@fsu.edu 16 17 Copyright 2002 Peter Beerli and Joseph Felsenstein, Seattle WA 18 Copyright 2003 Peter Beerli, Tallahassee FL 19 20 Permission is hereby granted, free of charge, to any person obtaining 21 a copy of this software and associated documentation files (the 22 "Software"), to deal in the Software without restriction, including 23 without limitation the rights to use, copy, modify, merge, publish, 24 distribute, sublicense, and/or sell copies of the Software, and to 25 permit persons to whom the Software is furnished to do so, subject 26 to the following conditions: 27 28 The above copyright notice and this permission notice shall be 29 included in all copies or substantial portions of the Software. 30 31 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 32 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 33 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 34 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR 35 ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF 36 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 37 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 38 39 $Id: migration.h 2165 2013-08-24 16:33:45Z beerli $ 40 *---------------------------------------------------------------- 41 */ 42 /** 43 * Definitions used for the program migrate 44 * this file holds structures and typedefs 45 * and defines 46 */ 47 48 49 /* Altivec support [NOT WORKING YET]*/ 50 #ifdef ALTIVEC 51 #include "altivec.h" 52 #endif /*ALTIVEC*/ 53 54 #ifdef ZNZ 55 #include "znzlib.h" 56 #endif 57 58 typedef char allele_type[DEFAULT_ALLELENMLENGTH]; //!< allele type string 59 60 typedef long twin_fmt[2]; 61 typedef long quad_fmt[4]; 62 typedef MYREAL pair[2]; 63 typedef long longpair[2]; 64 typedef float tetra[9]; 65 typedef float duo[2]; 66 #ifdef ALTIVEC 67 typedef FloatVec sitelike; 68 #else 69 #ifdef VARMUT 70 // best model for condlike? chunks of numcategs * (condlike_site * numsites) 71 // accessing a single 1D array? 72 // 73 typedef MYREAL * sitelike; //this allows to tread each site as a different datatype with its own mutation model 74 #else 75 #ifdef GAP 76 typedef MYREAL sitelike[5]; 77 #else 78 typedef MYREAL sitelike[4]; 79 #endif /*GAP*/ 80 #endif /*VARMUT*/ 81 #endif /*ALTIVEC*/ 82 83 typedef sitelike *ratelike; 84 typedef ratelike *phenotype; 85 typedef MYREAL contribarr[MAXCATEGS]; //!< \todo conditional likelihoods should change to 1-D vector 86 typedef short seqval[MAXCATEGS]; 87 88 typedef char allele_fmt[LINESIZE]; 89 //typedef char allele_fmt[DEFAULT_ALLELENMLENGTH]; 90 91 /// tipdate format 92 typedef struct tipdate_fmt 93 { 94 MYREAL date; 95 long id; 96 char * name; 97 } 98 tipdate_fmt; 99 100 /// records when the automatic stop of burnin happened 101 typedef struct _burnin_record_fmt 102 { 103 long locus; 104 long replicate; 105 long stopstep; 106 MYREAL ess; 107 MYREAL accept; 108 MYREAL variance; 109 MYREAL oldvariance; 110 long worker; 111 } burnin_record_fmt; 112 113 114 #ifdef BEAGLE 115 typedef struct _beagle_fmt { 116 boolean instance; 117 int *instance_handle; // indicator to the workorder for the GPU/CPU likelihood calculation 118 int numallocpartials; 119 int numpartials; 120 double *partials; // conditional likelihood for each site 121 int numallocallyweights; // pattern weighting 122 int *allyweights;// pattern weighting 123 double *weights; // these are the probabilities of the site rates (should move out to mutationmodels) 124 125 int numallocbranches; 126 int numbranches; 127 double *branch_lengths; // branchlengths for (for each edge leading to a child) 128 int *branch_indices; // indicator for transition matices 129 int numallocoperations; 130 int numoperations; 131 int *operations; // holds the operation for each parent,leftchild,lefttransitionmatrix,rightchild,rightt.. 132 133 int *scalingfactorsindices; 134 int scalingfactorscount; 135 136 } beagle_fmt; 137 #endif /*BEAGLE*/ 138 139 140 /// likelihood ratio test data structure 141 typedef struct lr_data_fmt 142 { 143 short type; //!< type of test \todo check whether this is obsolete 144 long elem; //!< number of elements in value1 and value2 145 char *value1; //!< values to test using a mix of numbers and m,c, etc 146 char *value2; //!< not used, reserved fo pairwise tests 147 char *connect;//!< custom migration matrix for value1 [no user interaction] 148 } 149 lr_data_fmt; 150 151 /// likelihood ratio test structure 152 typedef struct lratio_fmt 153 { 154 long alloccounter; //!< number of elements allocated in data 155 long counter; //!< number of elements in data 156 lr_data_fmt *data; //!< contains the LRT data 157 } 158 lratio_fmt; 159 160 /// temporary structure to hold precomputed values for the conditional likelihood calculation 161 typedef struct valrec 162 { 163 MYREAL rat, ratxi, ratxv, zz, z1, y1, ww1, zz1, ww2, zz2, z1zz, z1yy, xiz1, 164 xiy1xv, ww1zz1, vv1zz1, ww2zz2, vv2zz2; 165 } 166 valrec; 167 168 /// 2-D array to pointers to valrec structures 169 typedef valrec ***tbl_fmt; 170 171 /// union to hold the conditional likelihood data in tree nodes 172 typedef union xarray_fmt 173 { 174 MYREAL *a; //!< used for allelic types 175 phenotype s;//!< used for sequence data types 176 } 177 xarray_fmt; 178 179 #ifdef UEP 180 /// union to hold unique event polymorphism conditional likelihood data 181 typedef union ueparray_fmt 182 { 183 MYREAL *a; //!< holds probabilities for many state data \todo check the usage of UEP union 184 pair *s; //!< holds probabilities for 0/1 data 185 } 186 ueparray_fmt; 187 #endif 188 #ifdef VARMUT 189 typedef struct varmut_fmt 190 { 191 //sequence model 192 MYREAL freqa, freqt, freqg, freqc; //!< base frequencies 193 MYREAL freqr, freqy; //!< frequency of purins 194 MYREAL freqar, freqcy, freqgr, freqty; //!< frequency of transitions 195 MYREAL aa, bb; // 196 long oldsite; //!< number of site patterns (raw, with SNP adjustment 197 long endsite; //!< number of site patterns 198 MYREAL xi; //!< ratio of transitions 199 MYREAL xv; //!< ratio of transversions 200 MYREAL ttratio; //!< ration between transitions and transversions 201 MYREAL fracchange; // 202 //msat model 203 204 //allele model 205 206 } varmut_fmt; 207 #endif /*VARMUT*/ 208 209 210 /// sequence data structure 211 typedef struct seqmodel_fmt 212 { 213 // MYREAL freqa, freqt, freqg, freqc; //!< base frequencies 214 // MYREAL freqr, freqy; //!< frequency of purins 215 // MYREAL freqar, freqcy, freqgr, freqty; //!< frequency of transitions 216 // BASEFREQLENGTH is defined in definitions.h 217 MYREAL basefrequencies[BASEFREQLENGTH]; //!< base frequencies (a,c,g,t), gap, purins (r,y), transitions (ar,cy,gr,ty) 218 // rate matrix that should be in basefrequencies [g=nucleotide g, G=gap] 219 // pa pc pg pt pG 220 // rac rag rat raG rcg rct rcG rgt rgG rtG 221 // model vector (Huelsenbeck and Alfaro) augmented to allow the inclusion of gaps but using an additional class 0 --> set to zero. 222 // abcdefghij 223 // no gaps: 224 // minimal: 1110110100 225 // JC/F81: 1110110100 226 // F84: 1110110100 227 // maximal: 1230450600 228 // full model: 229 // minimal: 1111111111 230 // JC/F81: 1110110100 231 // maximal: 123456789a 232 // 233 MYREAL aa, bb; // 234 long oldsite; //!< number of site patterns (raw, with SNP adjustment 235 long endsite; //!< number of site patterns 236 MYREAL xi; //!< ratio of transitions 237 MYREAL xv; //!< ratio of transversions 238 MYREAL ttratio; //!< ration between transitions and transversions 239 MYREAL fracchange; // 240 long *sites; // 241 long *alias; // 242 long *ally; // 243 long *category; // 244 short *weight; // 245 long weightsum; // 246 long *aliasweight; // 247 long *location; // 248 long addon; // 249 boolean *links; //!< information about link status of locus 250 long *savealiasweight; //!< keeps a savecopy of the aliasweights 251 boolean done; 252 } 253 seqmodel_fmt; 254 255 /// defines the data structure read from infile 256 typedef struct _data 257 { 258 FILE *infile; //!< data file 259 FILE *utreefile; //!< user tree input file 260 FILE *weightfile; //!< site weighting file 261 FILE *catfile; //!< site categories file 262 FILE *sumfile; //!< intermediate output summary file 263 FILE *distfile; //!< distance file among individuals, instead of treefile 264 FILE *geofile; //!< geographic distance among sampling locations 265 FILE *divfile; //!< fixed divergence time among populations 266 FILE *datefile; //!< sample date for each individual 267 #ifdef UEP 268 FILE *uepfile; //!< unique event polymorphism file 269 int **uep; //!< uep data holder 270 long uepsites; //!< number of UEP sites 271 MYREAL uepfreq0; //!< base frequency for 0 272 MYREAL uepfreq1; //!< base frequency for 1 273 #endif 274 char *****yy; //!< data holder 275 MYREAL *geo; //!< geographic distance data holder 276 MYREAL *lgeo; //!< log of geo 277 MYREAL **ogeo; //!< original of geo as read from the geofile 278 char ***allele; // 279 long *maxalleles; //!< maximal number of alleles in dataset 280 boolean *skiploci; //!< number of loci to skip -- loci with no data 281 char **popnames; //!< array of the population names 282 long **numind; //!< array of the number of individuals 283 long **numalleles; // 284 tipdate_fmt ***sampledates; 285 MYREAL maxsampledate; 286 char ****indnames; //!< array of individual names per population 287 long numpop; //!< number of populations 288 long loci; //!< number of loci 289 long *position; //!< position on chromosome [what to do with mutiple chromosomes?] 290 seqmodel_fmt **seq; //!< sequence data model 291 MYREAL freq; // 292 MYREAL freqlast; // 293 char dlm; //!< delimiter for allelic data 294 boolean hasghost; //!< whether there is a population with no data or not \todo check if hasghost is still needed 295 long ***shuffled; //!< holds the random subset of individuals 296 boolean oneliner; //!< if a () is used in the sites lines it is assumed that all loci are on one line 297 char *datatype; 298 char *locitypes; //!< holds all datatypes for each locus, options->datatype is the default 299 char **locusname; 300 long *subloci; //!< holds the number of strictly linked "loci" that make up the locus for migrate 301 long numdatatypealloc; 302 long numsublocialloc; 303 long numlocitypesalloc; 304 long *numrepeatnumbers; 305 long **repeatnumbers; 306 long *repeatlength; 307 boolean has_repeats; 308 MYREAL *locusweight; //invariant loci treatment 309 } 310 data_fmt; 311 312 /// holds histogram and associated statistics for each locus 313 typedef struct _bayeshistogram 314 { 315 long binsum; //sum of all elements of all parameters and bins 316 long *bins; // number of bins for each parameter 317 float *results; // contains histogram, size is bins*numparam 318 char *set95; // holds in/out HPD set for 95% probability 319 char *set50; // holds in/out HPD set for 50% probability [this is attached to cred95] 320 // on a per parameter basis 321 // structure has a data storage vectors and the following are all pointers into it 322 long numparam; // number of parameters 323 MYREAL *datastore; // data storage, size is numparam*10 324 // pointers into data storage 325 MYREAL *minima; // contains minimal values for each parameter 326 MYREAL *maxima; // contains maximal values for each parameter 327 MYREAL *adjmaxima;// holds maxima values used in histogram [are smaller than maxima] 328 MYREAL *cred50l; // holds 50%-credibility margins (<all lower values>, 329 MYREAL *cred50u; //<all high values>) 330 MYREAL *cred95l; // holds 95%-credibility margins (<all lower values>) 331 MYREAL *cred95u; //<all high values>) 332 MYREAL *modes; // holds 95%-credibility margins (<all lower values>, <all high values>) 333 MYREAL *medians; 334 MYREAL *means; 335 MYREAL *stds; 336 MYREAL **covariance; // holds the covariance structure per locus 337 long n; //holds the number of samples taken, this seem to be missing for the reading from mdimfile, but available otherwise 338 } bayeshistogram_fmt; 339 340 /// holds data for Bayesian approach 341 typedef struct _bayes 342 { 343 long count; 344 long numpop2; 345 longpair *map; //mapping of parameter using custom migration matrix 346 MYREAL *datastore; //data store for 347 // pointers into datastore 348 MYREAL *priormean; // holds the means of priors 349 MYREAL *priorstd; // holds the standard deviation of priors [not always filled] 350 MYREAL *prioralpha; // holds alpha of prior distribution [not always filled] 351 MYREAL *delta; // change scale [how big can this be 352 MYREAL *minparam;// minimal allowed parameter value 353 MYREAL *maxparam; // maximal allowed parameter value 354 MYREAL *meanparam; // mean values for expprior 355 MYREAL *alphaparam; // alpha value for gamma 356 MYREAL *betaparam; // beta value for gamma 357 // record all changes of parameters 358 MYREAL *params; // save for parameter vectors 359 long allocparams; //number of allocated parameter vectors 360 long numparams; //number of saved parameter vectors 361 long paramnum; // which param one is working with 362 MYREAL starttime; // start time that was used to change the lineage on tree 363 MYREAL stoptime; // stop time that was used to change the lineage on the tree 364 MYREAL oldval; //holds prob(G|param) 365 long *datastore2; // datastore for long integer variables 366 // pointers into datastore2 367 long *accept; //holds numbers of accepted changes for each parameter AND the tree 368 long *trials; //holds how many time was tried to change a parameter 369 // holds histogram for each locus and summary statistics 370 bayeshistogram_fmt *histogram; 371 MYREAL *deltahist; // delta for all histograms per parameters [important for multilocus case 372 char * custm2; //pointer to the world->options->custm2; 373 boolean mu; 374 long mdiminterval ; // interval at which the whole parameter list is printed to file 375 int prettyhist; 376 MYREAL *scaling_factors; 377 MYREAL *histtotal; 378 MYREAL maxmaxvala; 379 // speed improvement to beat the slow malloc for large string allocation 380 // measures size of first string and uses that currently size * 2; 381 boolean has_linesize; 382 long linesize; 383 long progresslinesize; 384 long *mdimfilecount; 385 MYREAL *priors; 386 } 387 bayes_fmt; 388 389 /// storage for bayesian prior related parameters 390 typedef struct prior_fmt 391 { 392 MYREAL min; 393 MYREAL mean; 394 MYREAL std; 395 MYREAL max; 396 MYREAL delta; 397 MYREAL alpha; 398 long bins; 399 } prior_fmt; 400 401 402 /// data storage for all options 403 typedef struct _option 404 { 405 char **buffer; //pointer to buffer when MPI-worker otherwise NULL 406 FILE *parmfile; 407 FILE *seedfile; 408 FILE *logfile; 409 FILE *mixfile; 410 /*general options */ 411 // 412 // name length 413 long nmlength; /* length of individual names */ 414 long popnmlength; /* length of population names */ 415 long allelenmlength; /* length of allele names */ 416 // 417 // custom migration matrix setup 418 char *custm; /* custom migration matrix: theta are first */ 419 char *custm2; /*custom migration matrix: theta on diagonal */ 420 //long symn; 421 //long sym2n; 422 //long zeron; 423 //long constn; 424 //long tmn; 425 //long mmn; 426 // long *zeroparam; 427 // long *constparam; 428 // twin_fmt *symparam; 429 // quad_fmt *sym2param; 430 // long *mmparam; 431 432 /*input/output options */ 433 int menu; 434 boolean progress; 435 boolean verbose; 436 boolean writelog; 437 boolean uep; 438 MYREAL ueprate; 439 MYREAL uepmu; 440 MYREAL uepnu; 441 MYREAL uepfreq0; 442 MYREAL uepfreq1; 443 boolean movingsteps; 444 MYREAL acceptfreq; 445 boolean printdata; 446 boolean printcov; 447 boolean usertree; 448 boolean usertreewithmig; 449 boolean randomtree; 450 short treeprint; 451 long treeinc; 452 boolean usem; 453 boolean recordedusem; 454 short migvar; 455 boolean plot; 456 boolean plotnow; 457 short plotmethod; 458 short plotvar; 459 short plotscale; 460 MYREAL plotrange[4]; 461 long plotintervals; 462 MYREAL *plotxvalues; 463 MYREAL *plotyvalues; 464 boolean simulation; 465 boolean mighist; 466 boolean mighist_all; 467 long mighist_counter; 468 long mighist_increment; 469 boolean skyline; 470 float eventbinsize; 471 boolean mixplot; 472 boolean dist; 473 boolean geo; 474 boolean div; 475 char *parmfilename; 476 char *infilename; 477 char *outfilename; 478 char *pdfoutfilename; 479 char *logfilename; 480 char *mathfilename; 481 char *treefilename; 482 char *utreefilename; 483 char *catfilename; 484 char *weightfilename; 485 char *sumfilename; 486 char *mighistfilename; 487 char *skylinefilename; 488 char *distfilename; 489 char *geofilename; 490 char *divfilename; 491 char *bootfilename; 492 char *seedfilename; 493 char *mixfilename; 494 #ifdef UEP 495 char *uepfilename; 496 #endif 497 char *datefilename; 498 char *bayesfilename; 499 char *bayesmdimfilename; 500 boolean mdimdelete; 501 int use_compressed; 502 boolean prioralone; 503 504 char title[LINESIZE + 1]; 505 lratio_fmt *lratio; 506 boolean aic; 507 boolean fast_aic; 508 MYREAL aicmod; 509 FILE *aicfile; 510 char * aicfilename; 511 char aictype[3]; 512 char fsttype; 513 boolean printfst; 514 short profile; 515 char profilemethod; 516 long df; 517 boolean qdprofile; 518 boolean printprofsummary; 519 boolean printprofile; 520 short profileparamtype; 521 522 /* data options */ 523 char datatype; 524 boolean include_unknown; 525 short migration_model; 526 char dlm; 527 long micro_stepnum; 528 long micro_threshold; 529 int msat_option; 530 pair msat_tuning; 531 MYREAL **steps; 532 boolean interleaved; 533 MYREAL seqerror; 534 MYREAL *ttratio; 535 boolean fastlike; 536 boolean freqsfrom; 537 long categs; 538 MYREAL *rate; 539 long rcategs; 540 MYREAL *rrate; 541 MYREAL *probcat; 542 long seqrate_gamma_num; /*number of loci for gamma values*/ 543 MYREAL *seqrate_gamma; /* for gamma values for many loci*/ 544 MYREAL probsum; 545 boolean gammarates; 546 boolean autocorr; 547 MYREAL lambda; 548 boolean weights; 549 MYREAL freqa; 550 MYREAL freqc; 551 MYREAL freqg; 552 MYREAL freqt; 553 /* random number options */ 554 short autoseed; 555 unsigned long inseed; 556 unsigned long saveseed; 557 /* mcmc options */ 558 boolean bayes_infer; 559 boolean integrated_like; 560 /* slice sampling*/ 561 boolean slice_sampling[PRIOR_SIZE]; 562 MYREAL *slice_sticksizes; 563 MYREAL updateratio; 564 int bayesprior[3]; 565 int bayespretty; 566 prior_fmt *bayespriortheta; 567 prior_fmt *bayespriorm; 568 prior_fmt *bayespriorrate; 569 bayes_fmt *bayes; 570 boolean has_bayesfile; 571 boolean has_bayesmdimfile; 572 long bayesmdiminterval; 573 short thetaguess; 574 short migrguess; 575 boolean gamma; 576 boolean murates; 577 boolean murates_fromdata; 578 boolean bayesmurates; 579 MYREAL alphavalue; 580 MYREAL *mu_rates; 581 long muloci; 582 long *newpops; 583 long newpops_numalloc; 584 long newpops_numpop; 585 // boolean poprelabeled; 586 MYREAL *inheritance_scalars; 587 long inheritance_scalars_numalloc; 588 boolean has_datefile; 589 #ifdef LONGSUM 590 591 boolean fluctuate; 592 MYREAL *flucrates; 593 #endif 594 long burn_in; 595 char burnin_autostop; 596 short heating; 597 MYREAL *thetag; 598 long numthetag; 599 MYREAL *mg; 600 long nummg; 601 long schains; 602 long sincrement; 603 long ssteps; 604 long lchains; 605 long lincrement; 606 long lsteps; 607 MYREAL lcepsilon; 608 boolean gelman; 609 boolean gelmanpairs; 610 long pluschain; //how many chains are allowed , currently HARDCODED 611 boolean replicate; 612 long replicatenum; 613 long gridpoints; 614 long numpop; 615 MYREAL heat[50]; 616 boolean adaptiveheat; 617 long heating_interval; 618 long heated_chains; 619 /* save genealogy summary options */ 620 boolean readsum; 621 boolean writesum; 622 /* threading over loci */ 623 int cpu; 624 // 625 // fatal attraction of zero resistance 626 MYREAL minmigsumstat; 627 // 628 MYREAL *mutationrate_year; // keeps mutation rate per year for each locus 629 long mutationrate_year_numalloc; // 630 MYREAL generation_year; 631 long randomsubset; 632 unsigned long randomsubsetseed; 633 #ifdef SEASON 634 MYREAL *timings; // vector of time points (options) 635 long numtimings; // number of elements in the timings vector 636 MYREAL **timemodifiers; // vector of parameter modifier vectors 637 #endif 638 boolean heatedswap_off; 639 boolean has_autotune; 640 MYREAL autotune; 641 long totalsites; //for snp data 642 MYREAL *wattersons; 643 long *segregs; 644 boolean onlyvariable; // analyze only variable sequence loci 645 MYREAL locusweight; // use a weight to calculate the invariant locus 646 boolean has_variableandone; // use only one invariant locus and reweight 647 long firstinvariant; // locus number of first invariant locus to reweight using locusweight 648 } 649 option_fmt; 650 651 //defines the mutations models used in the extended "Ican do everything mutation model string thing" 652 typedef struct _mutationmodel { 653 long numpatterns; // number of site patterns 654 long numsites; 655 long numstates; // number of states in model: DNA=4, DNA+gap=5, msat>2 656 double *basefreqs; 657 long numsiterates; 658 double *siterates; 659 double *siteprobs; 660 char datatype; //specifices the model 661 double *eigenvectormatrix; 662 double *inverseeigenvectormatrix; 663 double *eigenvalues; 664 } mutationmodel_fmt; 665 666 /* used in the tree structure */ 667 /// defines a node in the tree structure 668 typedef struct _node 669 { 670 struct _node *next, *back; 671 boolean tip; 672 char type; 673 long number; 674 long pop; 675 long actualpop; 676 long id; 677 long bid; 678 boolean visited; 679 xarray_fmt x; 680 #ifdef UEP 681 682 int *uep; 683 ueparray_fmt ux; 684 MYREAL uepscale; 685 #endif 686 MYREAL *s; 687 MYREAL *scale; 688 char *nayme; 689 boolean top; 690 boolean dirty; 691 MYREAL v; 692 MYREAL tyme; 693 MYREAL length; 694 } 695 node; 696 697 /// defines a node in the the timelist structure with reference to the tree node 698 typedef struct vtlist 699 { 700 node *eventnode; /* node with age=tyme */ 701 MYREAL age; /* tyme from top nodelet */ 702 MYREAL interval; /* interval t[i+1].age - t[i].age */ 703 long *lineages; 704 long from; 705 long to; 706 /* long pop; */ 707 long slice; 708 boolean visited; 709 } 710 vtlist; 711 712 /// holds the tree, access should be done only through root 713 typedef struct tree_fmt 714 { 715 node **nodep; 716 node *root; 717 long pop; 718 long tips; 719 } 720 tree_fmt; 721 722 /// holds timeslist that is used in the MCMC run 723 typedef struct timelist_fmt 724 { 725 long numpop; 726 long copies; 727 long allocT; 728 long T; 729 long oldT; 730 vtlist *tl; 731 long *lineages; 732 } 733 timelist_fmt; 734 735 /// 736 /// holds minimal statistic for the longsum version, this allows to calculate changes in size etc through time 737 /// holds all tree information: ordered: list of lineages, frompop, topop, eventtime, eventtype 738 #ifdef LONGSUM 739 typedef struct longsum_fmt 740 { 741 long *lineages; 742 long *lineages2; 743 long fromto; //indicator into migration matrix mm2m(from,to) 744 long to; 745 MYREAL eventtime; 746 MYREAL interval; 747 char eventtype; 748 } 749 longsum_fmt; 750 #endif /*LONGSUM*/ 751 752 /// 753 /// time archive, holds minimal statistic for a single tree 754 /// is a compressed version of the longsum_fmt with many things precalculated 755 typedef struct tarchive_fmt 756 { 757 long copies; 758 MYREAL lcopies; 759 MYREAL *data; //hold content 760 #ifdef ALTIVEC 761 762 FloatVec *vdata; 763 #endif 764 765 MYREAL *point; // points into data 766 MYREAL *wait; // points into data 767 MYREAL *kt; //points to wait 768 MYREAL *km; //points to wait 769 MYREAL *p; // points to point 770 MYREAL *mindex; // points to point 771 #ifdef LONGSUM 772 773 longsum_fmt *longsum; // holds all tree information: ordered: list of lineages, frompop, topop, eventtime, eventtype 774 long longsumlen; // holds list of groups in longsum. 775 #endif /*LONGSUM*/ 776 777 } 778 tarchive_fmt; 779 780 /// 781 /// time archive holds all statistics for all trees 782 typedef struct timearchive_fmt 783 { 784 long allocT; 785 long T; 786 long numpop; 787 long sumtips; 788 MYREAL param_like; 789 MYREAL thb; 790 MYREAL alpha; 791 #ifdef ALTIVEC 792 // FloatVec *lcopiesvec; 793 // FloatVec *data; 794 #endif /*ALTIVEC*/ 795 796 tarchive_fmt *tl; // holds each tree (summary) 797 MYREAL *parameters; // holds all the parameter data 798 MYREAL *param; // pointer into parameters 799 MYREAL *param0; //pointer into parameters 800 MYREAL *lparam0; //pointer into parameters 801 MYREAL *likelihood; //pointer into parameters 802 long trials; 803 MYREAL normd; 804 } 805 timearchive_fmt; 806 807 808 /// holds pltting parameters 809 typedef struct _plotmax 810 { 811 MYREAL x1; 812 MYREAL y1; 813 MYREAL l1; 814 MYREAL x2; 815 MYREAL y2; 816 MYREAL l2; 817 } 818 plotmax_fmt; 819 820 /// holds quantile information 821 typedef struct _quantile_fmt 822 { 823 char *name; 824 MYREAL *param; 825 } 826 quantile_fmt; 827 828 /// format for migration events 829 /// holds time, from, to, and sumtips 830 ///typedef MYREAL migevent_fmt[4]; 831 typedef struct _migevent_fmt 832 { 833 float age; 834 int from; 835 int to; 836 long sumlines; 837 } migevent_fmt; 838 839 /// holds migration events for migration event time histogram 840 typedef struct _mighist_fmt 841 { 842 long allocsize; 843 long copies; 844 long weight; 845 long migeventsize; 846 migevent_fmt *migevents; 847 } 848 mighist_fmt; 849 850 /// holds all migration histograms for all loci 851 typedef struct _mighistloci_fmt 852 { 853 // records migration events into mighist container 854 mighist_fmt *mighist; 855 long mighistnum; 856 long allocsize; 857 858 // records all events parameters per bin derived from timelist 859 // should make communication over MPI much smaller 860 duo **migeventbins; 861 long *migeventbinnum; 862 // records all expected parameters per bin derived from timelist 863 tetra **eventbins; 864 long *eventbinnum; 865 MYREAL eventbinsize; // used for eventbins and migeventbins 866 867 } 868 mighistloci_fmt; 869 870 /// \todo what does this format do? 871 typedef struct _histogram_fmt 872 { 873 long count; 874 MYREAL *time; 875 long *weight; 876 } 877 histogram_fmt; 878 879 /// reduced set of options that are used to run the MCMC chain 880 typedef struct _worldoption 881 { 882 boolean gamma; 883 MYREAL alphavalue; 884 boolean murates; 885 boolean murates_fromdata; 886 long muloci; 887 MYREAL *mu_rates; 888 MYREAL *lmu_rates; 889 MYREAL *meanmu; 890 MYREAL *inheritance_scalars; 891 boolean prioralone; 892 MYREAL *heat; 893 MYREAL *averageheat; 894 #ifdef LONGSUM 895 896 boolean fluctuate; 897 #endif 898 899 short migration_model; 900 char *custm; 901 char *custm2; 902 MYREAL *thetag; 903 MYREAL *mg; 904 long zeron; 905 long *zeroparam; 906 long constn; 907 long *constparam; 908 long symn; 909 long sym2n; 910 twin_fmt *symparam; 911 quad_fmt *sym2param; 912 long tmn; 913 long mmn; 914 long *mmparam; 915 boolean mixplot; 916 boolean progress; 917 boolean writelog; 918 FILE *logfile; 919 boolean plotnow; 920 boolean verbose; 921 boolean replicate; 922 boolean gelman; 923 boolean gelmanpairs; 924 MYREAL lcepsilon; 925 boolean simulation; 926 char datatype; 927 long lsteps; 928 long lincr; 929 MYREAL loglsteps; 930 long treeprint; 931 long treeinc; 932 long movingsteps; 933 MYREAL acceptfreq; 934 long rcategs; 935 long categs; 936 short heating; 937 long heated_chains; 938 long heating_interval; 939 long heating_count; 940 boolean adaptiveheat; 941 char profilemethod; 942 boolean printprofile; 943 boolean printprofsummary; 944 long profileparamtype; 945 long df; 946 long lchains; 947 long replicatenum; 948 long *micro_threshold; 949 long micro_stepnum; 950 pair msat_tuning; 951 MYREAL ***steps; 952 MYREAL *rrate; 953 MYREAL *rate; 954 MYREAL *probcat; 955 long pluschain; 956 boolean mighist; 957 boolean mighist_all; 958 long mighist_counter; 959 long mighist_increment; 960 boolean skyline; 961 float eventbinsize; 962 long burn_in; 963 char burnin_autostop; 964 boolean usem; 965 short migvar; 966 boolean plot; 967 long plotmethod; 968 long plotintervals; 969 MYREAL plotrange[4]; 970 short plotscale; 971 long plotvar; 972 MYREAL *plotxvalues; 973 MYREAL *plotyvalues; 974 lratio_fmt *lratio; 975 boolean aic; 976 boolean fast_aic; 977 MYREAL aicmod; 978 FILE *aicfile; 979 char aictype[3]; 980 MYREAL lambda; 981 FILE *mixfile; 982 #ifdef UEP 983 984 boolean uep; 985 MYREAL ueprate; 986 MYREAL uepmu; 987 MYREAL uepnu; 988 MYREAL uepfreq0; 989 MYREAL uepfreq1; 990 #endif 991 992 boolean fastlike; 993 boolean bayes_infer; 994 boolean slice_sampling[PRIOR_SIZE]; 995 MYREAL *slice_sticksizes; 996 MYREAL updateratio; 997 boolean has_bayesfile; 998 boolean has_bayesmdimfile; 999 long bayesmdiminterval; 1000 MYREAL minmigsumstat; 1001 boolean has_datefile; 1002 MYREAL *mutationrate_year; // keeps mutation rate per year for each locus 1003 long mutationrate_year_numalloc; // 1004 MYREAL generation_year; 1005 boolean treeinmemory; 1006 boolean heatedswap_off; 1007 boolean has_autotune; 1008 MYREAL autotune; 1009 boolean onlyvariable; // analyze only variable sequence loci 1010 MYREAL locusweight; // use a weight to calculate the invariant locus 1011 boolean has_variableandone; // use only one invariant locus and reweight 1012 long firstinvariant; // locus number of first invariant locus to reweight using locusweight 1013 1014 } 1015 worldoption_fmt; 1016 1017 /// defines the time range of a unique event polymorphism 1018 typedef struct _ueptime 1019 { 1020 long size; 1021 long *populations; // population the nodes are in 1022 MYREAL *ueptime; // first elementis bottom most time , last elem is 1023 // topmost time,in between are for each migration events 1024 } 1025 ueptime_fmt; 1026 1027 /// reduced set of parameters related to the data used to run MCMC 1028 typedef struct _worlddata 1029 { 1030 boolean *skiploci; 1031 MYREAL *geo; 1032 MYREAL *lgeo; 1033 long *maxalleles; 1034 seqmodel_fmt **seq; 1035 FILE *sumfile; 1036 MYREAL freq; 1037 MYREAL freqlast; 1038 tipdate_fmt ***sampledates; 1039 MYREAL maxsampledate; 1040 long **numind; 1041 #ifdef UEP 1042 1043 long uepsites; 1044 #endif 1045 MYREAL *locusweight; //invariant loci treatment 1046 } 1047 worlddata_fmt; 1048 1049 //typedef struct _proposal_fmt proposal_fmt; 1050 1051 typedef struct _convergence 1052 { 1053 MYREAL *gelmanmeanmaxR; 1054 MYREAL gelmanmeanRall; 1055 MYREAL gelmanmaxRall; 1056 MYREAL *chain_s; 1057 MYREAL *chain_means; 1058 long *chain_counts; 1059 } convergence_fmt; 1060 1061 #ifdef MPI 1062 typedef struct _mpirequest_fmt 1063 { 1064 char *tempstr; 1065 int sender; 1066 int tag; 1067 } mpirequest_fmt; 1068 #endif 1069 /// holds all the parameters etc to run the program 1070 1071 //typedef struct proposal_fmt; 1072 1073 typedef struct _world 1074 { 1075 /* generalities */ 1076 char *name; 1077 char *worldname; 1078 worldoption_fmt *options; 1079 worlddata_fmt *data; 1080 1081 #ifdef MPI 1082 int *who; 1083 int *mpistack; 1084 long mpistack_numalloc; 1085 mpirequest_fmt *mpistack_request; 1086 long mpistack_request_numalloc; 1087 long mpistack_requestnum; 1088 long mpistacknum; 1089 #ifdef SLOWNET 1090 1091 int *profilewho; 1092 #endif 1093 #endif /*MPI*/ 1094 long allocbufsize; 1095 char *buffer; // buffer for profiles, not needed before profiletables() 1096 long loci; 1097 long skipped; /*loci with no data */ 1098 long locus; /* the current locus, if single then set to 0 */ 1099 long thislocus; /* the real current locus, the whole locus scheme 1100 needs revision */ 1101 long numpop; 1102 long numpop2; 1103 long sumtips; 1104 /* migration parameter array starts and ends */ 1105 int *mstart; 1106 int *mend; 1107 1108 /* time archives, contains the data/results for summarizing */ 1109 timearchive_fmt **atl; 1110 1111 /* migration histogram reporter */ 1112 mighistloci_fmt *mighistloci; 1113 long mighistlocinum; 1114 1115 /*tree material */ 1116 node **nodep; 1117 node *root; 1118 MYREAL treelen; 1119 long unique_id; 1120 tbl_fmt tbl; 1121 contribarr *contribution; 1122 /* parameter */ 1123 MYREAL *param0; 1124 MYREAL *param00; 1125 MYREAL **fstparam; 1126 /*allows for 3 different rate changes at specific times*/ 1127 #ifdef LONGSUM 1128 MYREAL *flucrates; 1129 MYREAL *lflucrates; 1130 #endif 1131 /* mcmc related */ 1132 long *lineages; 1133 timelist_fmt *treetimes; 1134 MYREAL *mig0list; 1135 MYREAL **migproblist; 1136 long *design0list; 1137 1138 /*nr related */ 1139 MYREAL ***apg0; /* part-loglikelihoods of param0 */ 1140 MYREAL ***apg; /* part-loglikelihoods of param */ 1141 1142 /*heating scheme */ 1143 long actualinc; 1144 long increment; 1145 // MYREAL heatratio; 1146 boolean cold; /* true for cold chain, false for all others*/ 1147 MYREAL heat; 1148 MYREAL averageheat; 1149 MYREAL varheat; 1150 long heatid; 1151 struct _proposal_fmt *proposal; 1152 boolean has_proposal_details; 1153 MYREAL essminimum; 1154 MYREAL logprior; 1155 long treeswapcount; 1156 /* reporting time */ 1157 time_t starttime; 1158 long treesdone; 1159 long treestotal; 1160 boolean in_burnin; 1161 burnin_record_fmt *burnin_stops; 1162 long burnin_z; 1163 /* gelman R reporting/checking "gelman-convergence option" */ 1164 convergence_fmt *convergence; 1165 1166 long chains; 1167 boolean start; //I need this in estimateParameter() 1168 /* reporting */ 1169 MYREAL *likelihood; /* data likelihoods */ 1170 long alloclike; /* how many likelihoods can we store */ 1171 long numlike; /* how many likelihoods maximally are stored */ 1172 plotmax_fmt **plotmax; 1173 quantile_fmt *quantiles; 1174 boolean **percentile_failed; /*indicator whether the percentile calculation failed*/ 1175 boolean percentile_some_failed; /*some precentiles failed*/ 1176 MYREAL maxdatallike; /* the maximum log likelihood of a chain */ 1177 MYREAL allikemax; /* the maximum log likelihood the best tree */ 1178 boolean in_last_chain; /*for last CHAIN option in print-tree */ 1179 MYREAL param_like; /* current parameter likelihood [=current chain] */ 1180 MYREAL **chainlikes; /* parameter likelihood of replicate/last chains all loci */ 1181 long trials; 1182 MYREAL normd; 1183 long repkind; 1184 long rep; // defines single cahin estimators 1185 long replicate; //holdss the replicate we are in 1186 long repstop; 1187 long lsteps; 1188 MYREAL ***cov; 1189 long migration_counts; 1190 char ****plane; 1191 #ifdef UEP 1192 1193 MYREAL ueplikelihood; 1194 MYREAL **ueplike; // contains uep like per tree pop x ueps 1195 MYREAL ****ueplikestore; // constains all ueplike long-steps x pop x uep 1196 ueptime_fmt *ueptime; //contains ueptime per tree and mutation 1197 ueptime_fmt ***ueptimestore; //contains all ueptimes 1198 long ***ueprootstore; //stores uep status at root 1199 long *oldrootuep; 1200 long **uepanc; 1201 #endif 1202 long accept; 1203 MYREAL accept_freq; 1204 long swapped; 1205 long G; 1206 long bayesaccept; 1207 bayes_fmt *bayes; 1208 FILE *bayesfile; 1209 #ifdef ZNZ 1210 znzFile bayesmdimfile; 1211 #else 1212 FILE *bayesmdimfile; 1213 #endif 1214 FILE *outfile; 1215 FILE *pdfoutfile; 1216 FILE *treefile; 1217 FILE *mathfile; 1218 FILE *mighistfile; 1219 FILE *skylinefile; 1220 char **treespace; 1221 long *treespacenum; 1222 long *treespacealloc; 1223 MYREAL *besttreelike; 1224 // Bayes factor 1225 MYREAL *bf; 1226 MYREAL *bfscale; 1227 MYREAL *hmscale; 1228 MYREAL *amscale; 1229 MYREAL *hm; 1230 MYREAL *am; 1231 // autocorrelation and ESS 1232 MYREAL *autocorrelation; 1233 MYREAL *effective_sample; 1234 MYREAL *auto_archive; 1235 MYREAL *ess_archive; 1236 1237 node **node_collection; 1238 long node_collection_count; 1239 long node_collection_allocated; 1240 #ifdef MPI 1241 #ifdef PARALIO 1242 MPI_File mpi_bayesmdimfile; 1243 #endif 1244 #endif 1245 long *nummutationmodels;//per locus 1246 long *sublocistarts;//per locus 1247 long *maxnumpattern;//per locus 1248 mutationmodel_fmt *mutationmodels; 1249 // 1250 float page_width; 1251 float page_height; 1252 // 1253 #ifdef BEAGLE 1254 beagle_fmt *beagle; 1255 #endif 1256 #ifdef SEASON 1257 MYREAL *timings; // vector of time points 1258 long numtimings; // number of elements in the timings vector 1259 MYREAL **timemodifiers; // vector of parameter modifier vectors 1260 #endif 1261 char * warning; 1262 long warningsize; 1263 long warningallocsize; 1264 long maxreplicate; 1265 long *accept_archive; 1266 long *trials_archive; 1267 } 1268 world_fmt; 1269 1270 1271 1272 /// helper for the ML maximizer, nr originally came from newton-raphson 1273 typedef struct _nr_fmt 1274 { 1275 long partsize; /*number of part-variables, fixed per model */ 1276 MYREAL *parts; /* parts of the first and second derivatives */ 1277 MYREAL *d; /* first derivates */ 1278 MYREAL *od; 1279 MYREAL *dv; 1280 MYREAL *delta; 1281 MYREAL *gdelta; 1282 MYREAL *param; /* changed values of param */ 1283 MYREAL *lparam; /* changed values of log(param) */ 1284 MYREAL *values; /* profile values */ 1285 MYREAL *locilikes; 1286 1287 // MYREAL *saved; /* saved first derivates (used in MPI code) */ 1288 /* migration parameter array starts and ends */ 1289 int *mstart; //link to world->mstart 1290 int *mend; //link to world->mend 1291 MYREAL normd; 1292 long repkind; 1293 long repstart; 1294 long repstop; 1295 // long numg; 1296 timearchive_fmt **atl; 1297 world_fmt *world; 1298 MYREAL **dd; /* second derivatives */ 1299 MYREAL llike; /* parameter LOGlikelihood */ 1300 MYREAL lastllike; /* parameter LOGlikelihood */ 1301 MYREAL ollike; /* old " " */ 1302 MYREAL *datalike; /*P(D|G) */ 1303 MYREAL ***apg0; /* part-loglikelihoods of param0 */ 1304 MYREAL ***apg; /* part-loglikelihoods of param */ 1305 MYREAL *apg_max; /* maxvalue of apg */ 1306 long seqrate_gamma_num; /*number of loci for gamma values*/ 1307 MYREAL *seqrate_gamma; /* for gamma values for many loci*/ 1308 MYREAL *rate; /* rates for gamma rates */ 1309 MYREAL *probcat; /* probabilities for gamma rates */ 1310 long categs; /* #categories for gamma rates */ 1311 MYREAL alpha; 1312 long numpop; /* number of populations */ 1313 long numpop2; /* 2*numpop */ 1314 MYREAL *PGC; /* uncorrected llike */ 1315 MYREAL *oPGC; /* uncorrected ollike */ 1316 long copy_nr; /* number of genealogies */ 1317 boolean *skiploci; 1318 long profilenum; /* number of profile parameters */ 1319 long *profiles; /* which profile parameters */ 1320 long *indeks; /* which noprofile parameters */ 1321 } 1322 nr_fmt; 1323 1324 /// helper format for the maximizer and the profile calculator 1325 typedef struct helper_fmt 1326 { 1327 long locus; 1328 nr_fmt *nr; 1329 timearchive_fmt **atl; 1330 //MYREAL *param; 1331 long which; 1332 MYREAL weight; 1333 MYREAL ll; 1334 boolean multilocus; 1335 boolean boolgamma; 1336 long analystype; 1337 MYREAL *dv; 1338 MYREAL *xv; 1339 MYREAL *expxv; 1340 MYREAL sign; 1341 MYREAL lamda; 1342 } 1343 helper_fmt; 1344 1345 1346 /// migration table 1347 typedef struct _migr_table_fmt 1348 { 1349 long from; 1350 long to; 1351 MYREAL time; 1352 } 1353 migr_table_fmt; 1354 1355 /// 1356 /// during the treeupdates we need a scratchpad, proposal_fmt 1357 /// holds all temporary data 1358 typedef struct _proposal_fmt 1359 { 1360 world_fmt *world; 1361 char datatype; 1362 long sumtips; 1363 long numpop; 1364 long endsite; 1365 MYREAL fracchange; 1366 MYREAL *param0; 1367 node *root; 1368 short migration_model; 1369 long listsize; // holds the size of several nodelists 1370 boolean mig_removed; 1371 MYREAL rr; 1372 node **nodedata; // holds abovenodes and bordernodes 1373 node *origin; // pointers into nodedata 1374 node *target; // .... 1375 node *realtarget; 1376 node *tsister; 1377 node *realtsister; 1378 node *osister; 1379 node *realosister; 1380 node *ocousin; 1381 node *realocousin; 1382 node *oback; 1383 node *realoback; 1384 node **line_f; 1385 node **line_t; 1386 node *connect; 1387 MYREAL likelihood; 1388 MYREAL ueplikelihood; 1389 MYREAL **ueplike; 1390 MYREAL time; 1391 MYREAL v; 1392 MYREAL vs; 1393 xarray_fmt xt; 1394 xarray_fmt xf; 1395 MYREAL *mf; 1396 MYREAL *mt; 1397 //could be used to calculate p(gn|go): MYREAL nu; 1398 #ifdef UEP 1399 ueparray_fmt ut; 1400 ueparray_fmt uf; 1401 node *firstuep; 1402 MYREAL *umf; 1403 MYREAL *umt; 1404 #endif 1405 1406 node **aboveorigin; 1407 node **bordernodes; 1408 migr_table_fmt *migr_table; 1409 migr_table_fmt *migr_table2; 1410 long migr_table_counter; 1411 long migr_table_counter2; 1412 long old_migr_table_counter; 1413 long old_migr_table_counter2; 1414 long timeslice; 1415 MYREAL *mig0list; 1416 long *design0list; 1417 MYREAL treelen; 1418 1419 #ifdef BEAGLE 1420 long parentid; 1421 long leftid; 1422 long rightid; 1423 #endif 1424 } proposal_fmt; 1425 1426 1427 /// splines are not yet used in the profiles 1428 typedef struct spline_fmt 1429 { 1430 long ntab; 1431 long nwork; 1432 MYREAL *param; 1433 MYREAL *like; 1434 MYREAL *diff; 1435 MYREAL *diff2; 1436 long *constr; 1437 long *diagn; 1438 MYREAL *work; 1439 } 1440 spline_fmt; 1441 1442 /// ????????????????? 1443 typedef struct locusdata_fmt 1444 { 1445 long locus; 1446 world_fmt **universe; 1447 option_fmt *options; 1448 data_fmt *data; 1449 } 1450 locusdata_fmt; 1451 1452 1453 /// holds the minimal statistic for AIC for each model 1454 typedef struct _aic_fmt 1455 { 1456 MYREAL aic; 1457 char *pattern; 1458 long numparam; 1459 MYREAL mle; 1460 MYREAL lrt; 1461 MYREAL prob; 1462 MYREAL probcorr; 1463 } 1464 aic_fmt; 1465 1466 /// holds all models that are used for the AIC 1467 typedef struct _aic 1468 { 1469 aic_fmt * aicvec; 1470 long aicnum; 1471 MYREAL *param0; 1472 } 1473 aic_struct; 1474 1475 /// 1476 /// plot options 1477 typedef struct _plotfield_fmt 1478 { 1479 boolean print; 1480 char type; // 'a' = ascii plot (see mig-histogram.c) , 'p' = PDF plotting (see pretty.c) 1481 long xsize; // width printpositions in chars 1482 long ysize; // height printpositions in chars 1483 char xaxis[255]; //xaxis label 1484 char yaxis[255]; //yaxis label 1485 char yfaxis[255]; //frequency yaxis label 1486 char title[255]; 1487 float *yfreq; 1488 long *y; 1489 char **data; //the plotplane 1490 } 1491 plotfield_fmt; 1492 1493 /// 3 longs used for random number 1494 typedef long longer[3]; /* used in random.c */ 1495 1496 1497 extern long *seed; 1498 extern long *newseed; 1499 #ifdef SLOWNET 1500 extern MYREAL (*calc_like) (helper_fmt *, MYREAL *, MYREAL *); 1501 extern void (*calc_gradient) (nr_fmt *, helper_fmt *, MYREAL *); 1502 extern void (*setup_param0) (world_fmt *, nr_fmt *, long, long, long, long, 1503 long, boolean); 1504 1505 #endif 1506 //#ifdef MEMDEBUG 1507 //#include <memcheck.h> 1508 //#endif 1509 #include "sort.h" 1510 //#include "migrate_mpi.h" 1511 #endif 1512