/*------------------------------------------------------ Maximum likelihood estimation of migration rate and effective population size using a Metropolis-Hastings Monte Carlo algorithm ------------------------------------------------------- M C M C R O U T I N E S Markov Monte Carlo stuff: treechange, acceptance Peter Beerli 1996, Seattle beerli@scs.fsu.edu Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA Copyright 2003-2006 Peter Beerli, Tallahassee FL This software is distributed free of charge for non-commercial use and is copyrighted. Of course, we do not guarantee that the software works and are not responsible for any damage you may cause or have. $Id: mcmc1.c 733 2007-05-31 02:28:34Z beerli $ -------------------------------------------------------*/ /* \file mcmc1.c Tree changer and acceptance rejection scheme */ #include "migration.h" #include "sighandler.h" #include "tools.h" #include "random.h" #include "tree.h" #include "mcmc2.h" #ifdef UEP #include "uep.h" #endif #include "bayes.h" #ifdef DMALLOC_FUNC_CHECK #include #endif #define MIGRATION_AIR (boolean) 1 #define MIGRATION_IN_TREE (boolean) 0 #define NO_MIGR_NODES 0 #define WITH_MIGR_NODES 1 extern int myID; /* global for this file */ proposal_fmt * gproposal; /* prototypes ------------------------------------------- */ // metropolize over trees long tree_update (world_fmt * world, long g); /* private functions */ void new_localtimelist (timelist_fmt ** ntl, timelist_fmt * otl, long numpop); void new_proposal (proposal_fmt ** proposal, timelist_fmt * tl, world_fmt * world); void set_new_proposal (proposal_fmt ** proposal, timelist_fmt * tl, world_fmt * world); void chooseOrigin (proposal_fmt * proposal); void construct_localtimelist (timelist_fmt * timevector, proposal_fmt * proposal); void traverseAllNodes (node * theNode, node *** nodelist, long *node_elem, long *oldnode_elem, int include_migration); void add_locallineages (timelist_fmt * timevector, proposal_fmt * proposal); void chooseTarget (proposal_fmt * proposal, timelist_fmt * timevector, node ** bordernodes, long *bordernum); void findbordernodes (node * theNode, proposal_fmt * proposal, long pop, node *** bordernodes, long *allocsize, long *bordernum, vtlist ** tyme, long gte); void free_proposal (proposal_fmt * proposal); void free_timevector (timelist_fmt * timevector); long prune_timelist (timelist_fmt * oldtv, timelist_fmt * newtv, node ** ptr, long numpop); long remove_doublettes_old (timelist_fmt * timevector, node ** ptr); long xor (node ** ptrl1, node ** ptrl2); long rmigrcount (proposal_fmt * proposal); int migrate_old (proposal_fmt * proposal, node * up, long *old_migr_table_counter, boolean air); int migrate (proposal_fmt * proposal, node * up); int migrateb (proposal_fmt * proposal, node * up); int pre_population (proposal_fmt * proposal, vtlist * ltime, long gte, long *slider); boolean acceptlike (world_fmt * world, proposal_fmt * proposal, long g, timelist_fmt * tyme); MYREAL eventtime (proposal_fmt * proposal, long pop, vtlist * tentry, char *event); node *showsister (node * theNode); void count_migrations (node * p, long *count); long migration_from (long to, proposal_fmt * proposal); MYREAL prob_tree (world_fmt * world, timelist_fmt * tyme); /* Functions ++++++++++++++++++++++++++++++++++++++++++++++++*/ void set_tree_dirty (node * p) { switch (p->type) { case 'm': set_dirty (p); set_tree_dirty (p->next->back); break; case 't': break; case 'i': set_dirty (p); set_tree_dirty (p->next->back); set_tree_dirty (p->next->next->back); break; case 'r': set_dirty (p); set_tree_dirty (p->next->back); break; } } /*=======================================================*/ long tree_update (world_fmt * world, long g) { /*return 1 if tree was acceped 0 otherwise */ static long treefilepos; /* write position in the treefile */ #ifdef TESTING2 static boolean have_proposal=FALSE; #endif boolean coalesced; // boolean test; char event; long slider; long bordernum; long actualpop = -99, zz; MYREAL endtime, nexttime, age; #ifdef UEP boolean uepsuccess = FALSE; #endif proposal_fmt *proposal; /*scratchpad in which all involved nodes are recorded and help arrays, e.g. migration arrays, are stored */ timelist_fmt *timevector = NULL; /* local timelist */ vtlist *tentry = NULL; /*pointer into timeslice */ /* --------------------------------------------------------- initialize local timelist and construct residual timelist find the start node and snip of the branch and node below */ #ifdef __MWERKS__ eventloop (); #endif new_localtimelist (&timevector, &world->treetimes[0], world->numpop); #ifdef TESTING2 if(have_proposal) set_new_proposal (&proposal, &world->treetimes[0], world); else { new_proposal (&gproposal, &world->treetimes[0], world); set_new_proposal (&proposal, &world->treetimes[0], world); have_proposal = TRUE; } #else new_proposal (&proposal, &world->treetimes[0], world); #endif chooseOrigin (proposal); // if(world->options->bayes_infer) { world->bayes->starttime = -proposal->origin->tyme; // we start with negative value to // indicate that we do not know yet whether we will accept this move or not // if the move is accepted then the time will be negated (->positiv), we record // all start and stop times one could easily imagine that some parameter values // never gets accepted by the genealogy, and when we record this for bayesallfile we would // like to know this. see further doen at the acceptance for the proposals } // construct_localtimelist (timevector, proposal); tentry = &(*timevector).tl[0]; // if (proposal->origin->type == 't') //{ // age = proposal->origin->tyme; //} //else //{ age = proposal->origin->tyme; zz = 0; while (tentry->age <= age && zz < (*timevector).T) { tentry = &(*timevector).tl[zz]; zz++; } //} nexttime = tentry->age; if ((*timevector).T > 1) endtime = (*timevector).tl[(*timevector).T - 2].age; else endtime = 0.0; proposal->time = age; coalesced = FALSE; /*------------------------------------ main loop: sliding down the tree */ slider = 0; #ifdef TESTINGDATE // printf("Time Pop Type Lineages\n"); #endif while (nexttime <= endtime) { actualpop = (proposal->migr_table_counter > 0) ? proposal->migr_table[proposal->migr_table_counter - 1].from : proposal->origin->pop; proposal->time = age + eventtime (proposal, actualpop, tentry, &event); #ifdef TESTINGDATE // printf("%0.5f %3li %c %3li",proposal->time,actualpop,event,tentry->lineages[0]); //for(long jj = 1;jj < proposal->world->numpop; jj++) //printf(" %3li",tentry->lineages[jj]); //printf("\n"); #endif // FPRINTF(stderr,"age=%f event=%c\n",proposal->time,event); if (proposal->time < nexttime) { if (event == 'm') { if (!migrate (proposal, proposal->origin)) { // we end up here when the migration events exceed the upper limit free_proposal (proposal); free_timevector (timevector); //warning("Event was migration event but could not match a lineage to it\n"); return 0; } age = proposal->time; continue; } else { /*coalesce */ chooseTarget (proposal, timevector, proposal->bordernodes, &bordernum); if(bordernum == 0) { warning("chooseTarget failed at age %f\n", proposal->time); free_proposal (proposal); free_timevector (timevector); return 0; } pretendcoalesce1p (proposal); #ifdef UEP uepsuccess = is_success_pseudo_uep (proposal); #endif coalesced = TRUE; break; } } /*end if proposal->time < nextime */ age = nexttime; tentry = &(*timevector).tl[(tentry->slice) + 1]; /*next entry in timelist */ nexttime = tentry->age; } if (!coalesced) { if (!pre_population (proposal, (*timevector).tl, (*timevector).T - 1, &slider)) { free_proposal (proposal); free_timevector (timevector); return 0; } // fprintf(stderr,"%i> last time=%f, end time = %f (%c)\n", myID, age,(*timevector).tl[(*timevector).T-1].eventnode->tyme,(*timevector).tl[(*timevector).T-1].eventnode->type ); pretendcoalesce1p (proposal); #ifdef UEP uepsuccess = is_success_pseudo_uep (proposal); #endif } if ( #ifdef UEP ((!world->options->uep && !uepsuccess) || (world->options->uep && uepsuccess)) && #endif (acceptlike (world, proposal, g, timevector))) { if (proposal->time > world->root->tyme) { /*saveguard */ world->root->tyme += proposal->time; } coalesce1p (proposal); //#ifdef TESTINGDATE //traverse_adjust(world->root->next->back, 1.0); //#endif // record the time interval that was used for the lineage. if(world->options->bayes_infer) { world->bayes->starttime = -world->bayes->starttime; world->bayes->stoptime = proposal->time; } // #ifdef UEP world->likelihood[g] = treelikelihood (world); if (world->options->uep) { update_uep (world->root->next->back, world); check_uep_root (world->root->next->back, world); world->treelen = 0.0; calc_treelength (world->root->next->back, &world->treelen); world->ueplikelihood = ueplikelihood (world); world->likelihood[g] = /*world->ueplikelihood +*/ world->likelihood[g]; } #else world->likelihood[g] = treelikelihood (world); #endif /* create a new timelist */ construct_tymelist (world, &world->treetimes[0]); if (world->options->treeprint != NONE && world->cold) print_tree (world, g, &treefilepos); world->migration_counts = 0; /* report the number of migration on the tree */ count_migrations (world->root->next->back, &world->migration_counts); free_proposal (proposal); free_timevector (timevector); return 1; /* new tree accepted */ } else { // record the time interval that was used for the lineage. if(world->options->bayes_infer) { world->bayes->stoptime = -proposal->time; } // } free_proposal (proposal); free_timevector (timevector); return 0; /* not accepted */ } #ifdef SLATKIN_IMPORTANCE long slatkin_importance(world_mt *world, long g) { long slice=0; // create timelist with times drawn using the parameters new_localtimelist (&timevector, &world->treetimes[0], world->numpop); while (slice < sumlineages) { actualpop = (proposal->migr_table_counter > 0) ? proposal->migr_table[proposal->migr_table_counter - 1].from : proposal->origin->pop; age = (*timevector).tl[slice].age; (*timevector).tl[slice].age = age + eventtime (proposal, actualpop, tentry, &event); } // calculate which two datapoints to join give the distance } #endif /*=======================================================*/ /// /// allocate the timelist, contains times and lineages at that time /// lineages have a large storage device that is accessed from the /// tl[i].lineages. void new_localtimelist (timelist_fmt ** ntl, timelist_fmt * otl, long numpop) { // long i; (*ntl) = (timelist_fmt *) mycalloc (1, sizeof (timelist_fmt)); (*ntl)->tl = (vtlist *) mymalloc ((*otl).allocT * sizeof (vtlist)); (*ntl)->allocT = otl->allocT; (*ntl)->T = otl->T; // not needed with direct approach without sort memcpy ((*ntl)->tl, otl->tl, otl->allocT * sizeof (vtlist)); allocate_lineages (ntl,(*ntl)->tl, (*ntl)->allocT, (*ntl)->T, numpop); // not needed with direct approach without sort memcpy ((*ntl)->lineages, otl->lineages, (*ntl)->allocT * numpop * sizeof (long)); } /// /// reuse a global timelist to store values, this should save /// malloc/free calls void new_localtimelist_new (timelist_fmt ** ntl, timelist_fmt * otl, long numpop) { // UNFINISHED // long i; // (*ntl) = (timelist_fmt *) mycalloc (1, sizeof (timelist_fmt)); if((*ntl)->allocT < otl->allocT) { (*ntl)->tl = (vtlist *) myrealloc ((*ntl)->tl,(*otl).allocT * sizeof (vtlist)); // memset((*ntl)->tl+((*ntl)->allocT),0,sizeof(vtlist)*((*otl)->allocT-(*ntl)->allocT)); (*ntl)->allocT = otl->allocT; } (*ntl)->T = otl->T; memcpy ((*ntl)->tl, otl->tl, otl->allocT * sizeof (vtlist)); allocate_lineages (ntl,(*ntl)->tl, (*ntl)->allocT, (*ntl)->T, numpop); memcpy ((*ntl)->lineages, otl->lineages, (*ntl)->allocT * numpop * sizeof (long)); } void new_proposal (proposal_fmt ** proposal, timelist_fmt * tl, world_fmt * world) { #ifdef UEP long j; #endif long mal = world->data->maxalleles[world->locus]; //long listsize = ((*tl).allocT + (*tl).T + 5); long listsize = 2*(world->sumtips + 2); long sumtips = world->sumtips; // allocate the scratchpad (contains a shadow of the tree) (*proposal) = (proposal_fmt *) mycalloc (1, sizeof (proposal_fmt)); (*proposal)->likelihood = -HUGE; // pointers and values to outside structures (*proposal)->world = world; (*proposal)->datatype = world->options->datatype; (*proposal)->sumtips = world->sumtips; (*proposal)->numpop = world->numpop; (*proposal)->endsite = world->data->seq->endsite; (*proposal)->fracchange = world->data->seq->fracchange; (*proposal)->param0 = world->param0; (*proposal)->root = world->root; (*proposal)->migration_model = world->options->migration_model; // precalculated values (*proposal)->mig0list = world->mig0list; (*proposal)->design0list = world->design0list; (*proposal)->listsize = listsize; // nodes above the picked node + migration nodes (*proposal)->aboveorigin = (node **) mycalloc (listsize, sizeof (node *)); // node data holding vector for bordernodes, line_f, line_t (*proposal)->nodedata = (node **) mycalloc (3 * listsize, sizeof (node *)); // adjacent nodes of the pcked nodes (*proposal)->bordernodes = (*proposal)->nodedata; // line_f nodes (*proposal)->line_f = (*proposal)->bordernodes + listsize; (*proposal)->line_t = (*proposal)->line_f + listsize; #ifdef checkdebugstuff // nodes above the picked node (*proposal)->aboveorigin = (node **) mycalloc (listsize, sizeof (node *)); // node data holding vector for aboveorigin, bordernodes, line_f, line_t // (*proposal)->nodedata = // (node **) mycalloc (2 * listsize + 2 * sumtips, sizeof (node *)); // adjacent nodes of the pcked nodes (*proposal)->bordernodes = (node **) mycalloc (listsize, sizeof (node *)); // line_f nodes (*proposal)->line_f = (node **) mycalloc (listsize, sizeof (node *));; (*proposal)->line_t = (node **) mycalloc (listsize, sizeof (node *));; #endif // mf holds also mt array (*proposal)->mf = (MYREAL *) mycalloc (2* (*proposal)->endsite, sizeof (MYREAL)); (*proposal)->mt = (*proposal)->mf + (*proposal)->endsite; if (strchr (SEQUENCETYPES, (*proposal)->datatype)) { allocate_xseq(&(*proposal)->xf, world->data->seq->endsite, world->options->rcategs); allocate_xseq(&(*proposal)->xt, world->data->seq->endsite, world->options->rcategs); } else { (*proposal)->xf.a = (MYREAL *) mycalloc (1, sizeof (MYREAL) * mal); (*proposal)->xt.a = (MYREAL *) mycalloc (1, sizeof (MYREAL) * mal); } (*proposal)->old_migr_table_counter = 4 * sumtips /* 100 */ ; (*proposal)->old_migr_table_counter2 = 4 * sumtips /* 100 */ ; (*proposal)->migr_table = (migr_table_fmt *) mycalloc (1, sizeof (migr_table_fmt) * (*proposal)->old_migr_table_counter); (*proposal)->migr_table2 = (migr_table_fmt *) mycalloc ((*proposal)->old_migr_table_counter2, sizeof (migr_table_fmt)); (*proposal)->migr_table_counter = 0; (*proposal)->migr_table_counter2 = 0; #ifdef UEP if (world->options->uep) { (*proposal)->ueplike = (MYREAL **) mycalloc (world->data->uepsites, sizeof (MYREAL *)); (*proposal)->ueplike[0] = (MYREAL *) mycalloc (world->numpop * world->data->uepsites, sizeof (MYREAL)); for (j = 1; j < world->data->uepsites; ++j) (*proposal)->ueplike[j] = (*proposal)->ueplike[0] + j * world->numpop; (*proposal)->ut.s = (pair *) mycalloc (world->data->uepsites, sizeof (pair)); (*proposal)->uf.s = (pair *) mycalloc (world->data->uepsites, sizeof (pair)); (*proposal)->umt = (MYREAL *) mycalloc (world->data->uepsites, sizeof (MYREAL)); (*proposal)->umf = (MYREAL *) mycalloc (world->data->uepsites, sizeof (MYREAL)); } #endif } /// /// sets new proposal structure using template gproposal, thins function will replace new_proposal() void set_new_proposal (proposal_fmt ** proposal, timelist_fmt * tl, world_fmt * world) { long oldsize=0; long newsize =0; long sumtips ; #ifdef UEP long j; #endif // long mal = world->data->maxalleles[world->locus]; long listsize = ((*tl).allocT + (*tl).T + 5); // long sumtips = world->sumtips; // allocate the scratchpad (contains a shadow of the tree) (*proposal) = gproposal; sumtips = (*proposal)->sumtips; // pointers and values to outside structures (*proposal)->world = world; (*proposal)->datatype = world->options->datatype; (*proposal)->sumtips = world->sumtips; (*proposal)->numpop = world->numpop; (*proposal)->endsite = world->data->seq->endsite; (*proposal)->fracchange = world->data->seq->fracchange; (*proposal)->param0 = world->param0; (*proposal)->root = world->root; (*proposal)->migration_model = world->options->migration_model; // precalculated values (*proposal)->mig0list = world->mig0list; (*proposal)->design0list = world->design0list; if(listsize > (*proposal)->listsize) { newsize = 2 * listsize + 2 * sumtips; (*proposal)->nodedata = (node **) myrealloc ((*proposal)->nodedata, newsize * sizeof (node *)); oldsize = (2*((*proposal)->listsize)+2*sumtips); memset((*proposal)->nodedata+oldsize, 0, (newsize-oldsize) * sizeof (node *)); (*proposal)->aboveorigin = (*proposal)->nodedata; (*proposal)->bordernodes = (*proposal)->nodedata + listsize; (*proposal)->line_f = (*proposal)->nodedata + listsize + sumtips; (*proposal)->line_t = (*proposal)->nodedata + listsize + 2 * sumtips; } (*proposal)->listsize = listsize; // mf holds also mt array //(*proposal)->mf = (MYREAL *) mycalloc (2* (*proposal)->endsite, sizeof (MYREAL)); //(*proposal)->mt = (*proposal)->mf + (*proposal)->endsite; //if (strchr (SEQUENCETYPES, (*proposal)->datatype)) //{ // allocate_xseq(&(*proposal)->xf, world->data->seq->endsite, world->options->rcategs); // allocate_xseq(&(*proposal)->xt, world->data->seq->endsite, world->options->rcategs); //} //else //{ // (*proposal)->xf.a = (MYREAL *) mycalloc (1, sizeof (MYREAL) * mal); // (*proposal)->xt.a = (MYREAL *) mycalloc (1, sizeof (MYREAL) * mal); // } //(*proposal)->old_migr_table_counter = 4 * sumtips /* 100 */ ; // (*proposal)->old_migr_table_counter2 = 4 * sumtips /* 100 */ ; //(*proposal)->migr_table = // (migr_table_fmt *) mycalloc (1, // sizeof (migr_table_fmt) * // (*proposal)->old_migr_table_counter); // (*proposal)->migr_table2 = // (migr_table_fmt *) mycalloc ((*proposal)->old_migr_table_counter2, // sizeof (migr_table_fmt)); (*proposal)->migr_table_counter = 0; (*proposal)->migr_table_counter2 = 0; #ifdef UEP if (world->options->uep) { // (*proposal)->ueplike = // (MYREAL **) mycalloc (world->data->uepsites, sizeof (MYREAL *)); // (*proposal)->ueplike[0] = // (MYREAL *) mycalloc (world->numpop * world->data->uepsites, // sizeof (MYREAL)); // for (j = 1; j < world->data->uepsites; ++j) // (*proposal)->ueplike[j] = (*proposal)->ueplike[0] + j * world->numpop; (*proposal)->ut.s = (pair *) mycalloc (world->data->uepsites, sizeof (pair)); (*proposal)->uf.s = (pair *) mycalloc (world->data->uepsites, sizeof (pair)); (*proposal)->umt = (MYREAL *) mycalloc (world->data->uepsites, sizeof (MYREAL)); (*proposal)->umf = (MYREAL *) mycalloc (world->data->uepsites, sizeof (MYREAL)); } #endif } void jumblenodes (node ** s, long n) { node **temp; long i, rr, tn = n; temp = (node **) mycalloc (1, sizeof (node *) * n); memcpy (temp, s, sizeof (node *) * n); for (i = 0; i < n && tn > 0; i++) { s[i] = temp[rr = RANDINT (0, tn - 1)]; temp[rr] = temp[tn - 1]; tn--; } myfree(temp); } void traverse_check(node *theNode) { if (theNode == NULL) error("no node?????????"); if (theNode->type != 't') { if (theNode->next->back != NULL) { if(theNode->tyme < showtop(theNode->next->back)->tyme) { printf("Problem in traverse_check: id=%li time=%f type=%c time-up-next=%f type_up=%c\n", theNode->id, theNode->tyme, theNode->type, theNode->next->back->tyme, theNode->next->back->type); error("time conflict"); } traverse_check (theNode->next->back); } if (theNode->type != 'm' && theNode->next->next->back != NULL) { if(theNode->tyme < theNode->next->next->back->tyme) { printf("Problem in traverse_check: id=%li time=%f type=%c time-up-next-next=%f type_up=%c\n", theNode->id, theNode->tyme, theNode->type, theNode->next->next->back->tyme, theNode->next->back->type); error("time conflict"); } traverse_check (theNode->next->next->back); } } } /// /// pick a node at random from the available list of nodes, this code will even /// work with datasets where there are only two nodes. void chooseOrigin (proposal_fmt * proposal) { long elem = 0, oldelem = (proposal->sumtips * 2.); node *tmp=NULL, **goal; goal = (node **) mycalloc (oldelem, sizeof (node *)); //traverse_check(crawlback (proposal->root->next)); traverseAllNodes (crawlback (proposal->root->next), &goal, &elem, &oldelem, NO_MIGR_NODES); switch(elem) { case 0: error("problem with choosing a node for the tree change [choosOrigin()]"); break; case 1: case 2: tmp = goal[0]; break; default: tmp = goal[RANDINT (0, elem - 2)]; break; } myfree(goal); proposal->origin = tmp; if (proposal->origin != showtop (crawlback (proposal->root->next))) { proposal->oback = showtop (crawlback (proposal->origin)); proposal->osister = showsister (proposal->origin); if (proposal->oback != showtop (crawlback (proposal->root->next))) { proposal->ocousin = showsister (proposal->oback); } else { proposal->ocousin = NULL; } } if (proposal->origin == NULL) error ("Designation of origin for branch removal failed"); } void construct_localtimelist (timelist_fmt * timevector, proposal_fmt * proposal) { long z = 0; long oz = proposal->listsize; // long numpop = timevector->numpop = proposal->numpop; traverseAllNodes (crawlback (proposal->origin)->back, &proposal->aboveorigin, &z, &oz, WITH_MIGR_NODES); proposal->aboveorigin[z++] = proposal->oback; proposal->aboveorigin[z] = NULL; //z = remove_doublettes_old (timevector, proposal->aboveorigin); z = prune_timelist(&proposal->world->treetimes[0],timevector, proposal->aboveorigin, proposal->world->numpop); // not needed with direct approach without sort //qsort ((void *) (*timevector).tl, (*timevector).T, sizeof (vtlist), agecmp); //(*timevector).T -= z; if ((*timevector).tl[(*timevector).T - 1].eventnode->type != 'r') { error ("Root not at the end of timelist\n"); } timeslices (&timevector); add_locallineages (timevector, proposal); } /*---------------------------------------------------------------------------- finds all nodes in a tree starting at the root node and crawling up to the tips in a recursive fashion, writing nodeptrs in the nodelist vector the flag include_migration is 1 if we want to touch the migration nodes too, otherwise =0 -> jump over the migration nodes. for convenience we define the the macros NO_MIGR_NODES=0 and WITH_MIGR_NODES=1 in the treesetup.h file PB 1995 */ void traverseAllNodes (node * theNode, node *** nodelist, long *node_elem, long *oldnode_elem, int include_migration) { static long counter=0; long elem; counter++; if(theNode->type != 'r' && theNode->tyme > showtop(theNode->back)->tyme) { printf("function was called %li times\n",counter); error("traverseAllNodes() failed"); } if (include_migration == NO_MIGR_NODES) { // nodelist needs to be at least twice as long as sumtips // otherwise it will run out of reserved memory if (theNode->type != 't') { if (crawlback (theNode->next) != NULL) traverseAllNodes (crawlback (theNode->next), nodelist, node_elem, oldnode_elem, NO_MIGR_NODES); if (theNode->type != 'm' && crawlback (theNode->next->next) != NULL) traverseAllNodes (crawlback (theNode->next->next), nodelist, node_elem, oldnode_elem, NO_MIGR_NODES); (*nodelist)[(*node_elem)] = theNode; (*node_elem) += 1; if (theNode->type == 'm') { error ("Migration node encountered?! and died!"); } if(*node_elem > *oldnode_elem) error("die on the spot, this should not happen. Failed in traverseAllNodes without migration nodes"); } else { (*nodelist)[(*node_elem)] = theNode; (*node_elem) += 1; if(*node_elem > *oldnode_elem) error("die on the spot, this should not happen. Failed in traverseAllNodes without migration nodes at tip"); } } else { if (theNode->type != 't') { if (theNode->next->back != NULL) traverseAllNodes (theNode->next->back, nodelist, node_elem, oldnode_elem, WITH_MIGR_NODES); if (theNode->type != 'm' && theNode->next->next->back != NULL) traverseAllNodes (theNode->next->next->back, nodelist, node_elem, oldnode_elem, WITH_MIGR_NODES); if ((*node_elem) == (*oldnode_elem - 2)) { elem = *oldnode_elem = ((*oldnode_elem) + (*oldnode_elem) / 2); (*nodelist) = (node **) myrealloc ((*nodelist), sizeof (node *) * elem); memset ((*nodelist) + (*oldnode_elem), 0, sizeof (node *) * (elem - (*oldnode_elem))); *oldnode_elem = elem; } (*nodelist)[(*node_elem)++] = theNode; } else { if ((*node_elem) == (*oldnode_elem - 2)) { elem = *oldnode_elem = ((*oldnode_elem) + (*oldnode_elem) / 2); (*nodelist) = (node **) myrealloc ((*nodelist), sizeof (node *) * elem); memset ((*nodelist) + (*oldnode_elem), 0, sizeof (node *) * (elem - (*oldnode_elem))); *oldnode_elem = elem; } (*nodelist)[(*node_elem)++] = theNode; } } } /// /// copies elements from the old timelist (oldtv) /// into the new local timeslist (newtv) checking whether they /// are used after the removal of the residual tree that is above the origin (ptr) long prune_timelist (timelist_fmt * oldtv, timelist_fmt * newtv, node ** ptr, long numpop) { long i = 0; long j = 0; long slot = 0; long size = numpop * sizeof(long); // vtlist *otl = (*oldtv).tl; //vtlist *ntl = (*newtv).tl; long * tmpptr = NULL; node *thenode; //printf("oldtv.T=%li newtv.T=%li ", (*oldtv).T, (*newtv).T); /* assumes that there is an NULL element at the end */ slot=0; for (i = 0; i < (*oldtv).T; i++) { j = 0; thenode = (*oldtv).tl[i].eventnode; while ((thenode != ptr[j]) && (ptr[j] != NULL)) { // printf("%5li %5li old=%p ptr=%p\n",i, j ,(*oldtv).tl[i].eventnode,ptr[j]); j++; } // printf("*** %5li %5li old=%p ptr=%p\n",i, j ,(*oldtv).tl[i].eventnode,ptr[j]); if (ptr[j] == NULL) { tmpptr = (*newtv).tl[slot].lineages; memcpy(&(*newtv).tl[slot],&(*oldtv).tl[i],sizeof(vtlist)); (*newtv).tl[slot].lineages = tmpptr; memcpy((*newtv).tl[slot].lineages,(*oldtv).tl[i].lineages,size); //printf("in tree %li %li (%li)\n",z,i,otl[i].eventnode->number); slot++; } // else // { //printf("excised from timelist %li %li\n", i, otl[i].eventnode->number); // (*newtv).T -= 1; // } } (*newtv).T = slot; //printf("oldtv.T=%li newtv.T=%li last type=%c\n", (*oldtv).T, (*newtv).T, (*newtv).tl[(*newtv).T-1].eventnode->type); return (*oldtv).T - slot; } /* prepares to remove elements of a timelist by setting eventnode to NULL and age to Infinity, so a simple sort will push them over the edge */ long remove_doublettes_old(timelist_fmt * timevector, node ** ptr) { long i = 0, j = 0, slot = 0; /* assumes that there is an NULL element at the end */ for (i = 0; i < (*timevector).T; j = 0, i++) { while (((*timevector).tl[i].eventnode != ptr[j]) && (ptr[j] != NULL)) j++; if (ptr[j] != NULL) { slot++; (*timevector).tl[i].age = DBL_MAX; } } return slot; } void add_locallineages (timelist_fmt * timevector, proposal_fmt * proposal) { long pop = -99, numpop = proposal->world->numpop; if (timevector->T <= 0) error ("Help: timelist contains 0 elements"); for (pop = 0; pop < numpop; pop++) timevector->tl[timevector->T - 1].lineages[pop] = 0; if (timevector->T == 1) timevector->tl[0].lineages[proposal->osister->actualpop] += 1; else timevector->tl[timevector->T - 1].lineages[timevector->tl[timevector->T-2].from] += 1; add_partlineages (numpop, &timevector); } /* replaces nodepointers in list 1 with NULL if they are present in list 2 returns the first NULL slot in the array. */ long xor (node ** ptrl1, node ** ptrl2) { long i = 0, j = 0, slot = -1; /* assumes that there is an NULL element at the end */ for (i = 0; ptrl1[i] != NULL; j = 0, i++) { while ((ptrl1[i] != ptrl2[j]) && (ptrl2[j] != NULL)) j++; if (ptrl2[j] != NULL) { if (slot == -1) slot = i; ptrl1[i] = NULL; } } return slot; } /* migrate() fills the PROPOSAL->MIGRATION_TABLE in the tree or PROPOSAL->MIGRATION_TABLE2 when at the bottom of the tree PROPOSAL proposal-scratchpad UP node above (younger) MIGR_TABLE_COUNTER migration array counter, will increase by one during execution AIR if true standard execution, if false updating the last lineage in the residual tree. */ int migrate (proposal_fmt * proposal, node * up) { long i = proposal->migr_table_counter; migr_table_fmt *array = proposal->migr_table; if (i > MIGRATION_LIMIT * proposal->numpop) { return 0; } if (i > 0) array[i].to = array[i - 1].from; else array[i].to = up->pop; array[i].from = migration_from (array[i].to, proposal); //DEBUG // printf("i: %3li %li %li\n",i,proposal->migr_table[i].from,proposal->migr_table[i].to); array[i++].time = proposal->time; if (i >= proposal->old_migr_table_counter) { proposal->old_migr_table_counter += 10; proposal->migr_table = (migr_table_fmt *) myrealloc (proposal->migr_table, sizeof (migr_table_fmt) * (proposal->old_migr_table_counter)); } proposal->migr_table_counter = i; return 1; } int migrateb (proposal_fmt * proposal, node * up) { long i = proposal->migr_table_counter2; migr_table_fmt *array = proposal->migr_table2; if (i > MIGRATION_LIMIT * proposal->numpop) { return 0; } if (i > 0) array[i].to = array[i - 1].from; else array[i].to = up->pop; array[i].from = migration_from (array[i].to, proposal); // printf("t: %3li %li %li\n",i,proposal->migr_table2[i].from,proposal->migr_table2[i].to); array[i].time = proposal->time; i++; if (i >= proposal->old_migr_table_counter2) { proposal->old_migr_table_counter2 += 10; proposal->migr_table2 = (migr_table_fmt *) myrealloc (proposal->migr_table2, sizeof (migr_table_fmt) * (proposal->old_migr_table_counter2)); } proposal->migr_table_counter2 = i; return 1; } int migrate_old (proposal_fmt * proposal, node * up, long *old_migr_table_counter, boolean air) { migr_table_fmt *array; long i; if (air) { array = proposal->migr_table; i = proposal->migr_table_counter; } else { array = proposal->migr_table2; i = proposal->migr_table_counter2; } if (i > MIGRATION_LIMIT * proposal->numpop) { // FPRINTF (stdout, "migration limit reached\n"); return 0; } switch (proposal->migration_model) { case ISLAND: case ISLAND_VARTHETA: case MATRIX: case MATRIX_SAMETHETA: case MATRIX_ARBITRARY: if (i > 0) array[i].to = array[i - 1].from; else array[i].to = up->pop; array[i].from = migration_from (array[i].to, proposal); // printf("i: %3li %li %li\n",i,proposal->migr_table[i].from,proposal->migr_table[i].to); // printf("b: %3li %li %li\n",i,proposal->migr_table2[i].from,proposal->migr_table2[i].to); break; case CONTINUUM: case STEPSTONE: error ("not yet implemented\n"); break; default: break; } array[i++].time = proposal->time; if (i >= (*old_migr_table_counter)) { (*old_migr_table_counter) += 10; if (air) { proposal->migr_table = (migr_table_fmt *) myrealloc (proposal->migr_table, sizeof (migr_table_fmt) * (*old_migr_table_counter)); array = proposal->migr_table; } else { proposal->migr_table2 = (migr_table_fmt *) myrealloc (proposal->migr_table2, sizeof (migr_table_fmt) * (*old_migr_table_counter)); array = proposal->migr_table; } } if (air) { proposal->migr_table_counter = i; } else { proposal->migr_table_counter2 = i; } return 1; } /* migration_from() returns the FROM population when there was a migration TO population to migrate to PROPOSAL proposal-scratchpad */ long migration_from_old (long to, proposal_fmt * proposal) { long j, ii, msta, msto; MYREAL *geo = proposal->world->data->geo; MYREAL *r, rr = UNIF_RANDUM (); r = (MYREAL *) mycalloc (1, sizeof (MYREAL) * proposal->numpop); msta = mstart (to, proposal->numpop); msto = mend (to, proposal->numpop); r[0] = proposal->param0[msta] * geo[msta]; for (j = 1, ii = msta + 1; ii < msto; j++, ii++) { r[j] = r[j - 1] + geo[ii] * proposal->param0[ii]; } ii = 0; while (rr > r[ii] / r[j - 1]) { ii++; } myfree(r); if (ii < to) return ii; else return ++ii; } long migration_from (long to, proposal_fmt * proposal) { long ii = 0; MYREAL *r = proposal->world->migproblist[to]; MYREAL rr = UNIF_RANDUM (); while (rr > r[ii]) { ii++; } if (ii < to) return ii; else return ++ii; } void chooseTarget (proposal_fmt * proposal, timelist_fmt * timevector, node ** bordernodes, long *bordernum) { long actualpop = -99; node *rb = crawlback (proposal->root->next); *bordernum = 0; proposal->target = NULL; proposal->realtarget = NULL; if (proposal->migr_table_counter == 0) actualpop = proposal->origin->pop; else actualpop = proposal->migr_table[proposal->migr_table_counter - 1].from; //if (rb->tyme < proposal->time) //{ // error ("Wrong Time for action in chooseTarget()\n"); // proposal->target = NULL; // proposal->tsister = NULL; // proposal->realtsister = NULL; // proposal->realtarget = NULL; // return; //} findbordernodes (rb, proposal, actualpop, &bordernodes, &proposal->listsize, bordernum, &(*timevector).tl, (*timevector).T); if (*bordernum > 0) { proposal->target = bordernodes[RANDINT (0, (*bordernum) - 1)]; if (proposal->target != rb) { proposal->tsister = showsister (proposal->target); proposal->realtsister = crawlback (proposal->tsister)->back; } else proposal->tsister = NULL; proposal->realtarget = proposal->target; if (proposal->target->type == 'm') proposal->target = crawlback (showtop (proposal->target)->next); } else { proposal->target = NULL; proposal->tsister = NULL; proposal->realtsister = NULL; proposal->realtarget = NULL; } } void findbordernodes (node * theNode, proposal_fmt * proposal, long pop, node *** bordernodes, long *allocsize, long *bordernum, vtlist ** tyme, long gte) { node *tmp, *back; if (theNode == proposal->oback) { tmp = showtop (crawlback (proposal->osister)->back); back = showtop (proposal->oback->back); } else { tmp = showtop (theNode); back = showtop (theNode->back); } if (pop == tmp->pop && pop == back->actualpop && tmp->tyme < proposal->time && back->tyme > proposal->time) { if(*bordernum >= *allocsize) { warning("reallocation in finderbordernodes -- strange"); (*bordernodes) = (node **) myrealloc(*bordernodes, (*allocsize+50) * sizeof(node*)); memset(*bordernodes + *allocsize, 0, sizeof(node*) * 50); *allocsize = *allocsize + 50; } (*bordernodes)[(*bordernum)++] = tmp; return; } else { if (back->tyme < proposal->time) return; if (tmp->type != 't') { if (tmp->next->back != NULL) findbordernodes (tmp->next->back, proposal, pop, bordernodes, allocsize, bordernum, tyme, gte); if (tmp->type != 'm' && tmp->next->next->back != NULL) findbordernodes (tmp->next->next->back, proposal, pop, bordernodes, allocsize, bordernum, tyme, gte); } } } /* boolean same_pop(node * up, MYREAL tyme, long pop) { node *oldnn = showtop(up->back); node *nn = up; while (nn->tyme < tyme) { oldnn = nn; nn = showtop(nn->back); } if (oldnn->pop == pop && nn->actualpop == pop) return TRUE; else return FALSE; } */ /* ----------------------------------------------------------------------- simulates two lineages at once, if we are extending below the last node */ int pre_population (proposal_fmt * proposal, vtlist * ltime, long gte, long *slider) { boolean coalesced = FALSE; boolean choice = FALSE; long pop1 = -99, pop2 = -98; // long msta1=0, msto1=0; // long msta2=0, msto2=0; MYREAL ux; MYREAL rate = proposal->world->options->mu_rates[proposal->world->locus]; MYREAL invrate = 1./rate; MYREAL age1, denom, rr, r0, r1, horizon, mm, mm2; // MYREAL denom_invrate; if (gte > 0) proposal->realtarget = ltime[gte - 1].eventnode; else proposal->realtarget = ltime[0].eventnode->next->back; //?????gte if (proposal->realtarget == proposal->oback) { proposal->realtarget = crawlback (proposal->osister)->back; } if (proposal->realtarget->type == 'm') { proposal->target = crawlback (proposal->realtarget->next); if (proposal->target == proposal->oback) { proposal->target = proposal->osister; } } else { proposal->target = proposal->realtarget; } proposal->tsister = NULL; pop2 = proposal->realtarget->pop; pop1 = proposal->migr_table_counter > 0 ? proposal->migr_table[proposal->migr_table_counter - 1].from : proposal->origin->pop; age1 = MAX (proposal->realtarget->tyme, proposal->migr_table_counter > 0 ? proposal->migr_table[proposal->migr_table_counter - 1].time : proposal->origin->tyme); horizon = MAX (proposal->oback->tyme, age1); while (age1 < horizon) { mm = proposal->mig0list[pop1] * invrate; //mm = 0.0; //msta1 = mstart(pop1,proposal->numpop); //msto1 = mend(pop1,proposal->numpop); // for (i = msta1; i < msto1; i++) //{ //mm += proposal->param0[i]; //} if (pop1 == pop2) { denom = mm + (2. / (proposal->param0[pop1] * rate)); // denom_invrate = denom * invrate; ux = UNIF_RANDUM(); // proposal->time = age1 - LOG(ux) / denom_invrate; proposal->time = age1 - LOG(ux) / denom; age1 = proposal->time; if (age1 < horizon) { rr = UNIF_RANDUM (); r0 = (2. / (proposal->param0[pop1] * rate)) / denom; if (rr < r0) { return 1; } } } else { denom = mm; proposal->time = age1 - LOG (UNIF_RANDUM ()) / denom; age1 = proposal->time; } if (age1 < horizon) { if (!migrate (proposal, proposal->origin)) // &proposal->old_migr_table_counter, MIGRATION_AIR)) { return 0; } pop1 = proposal->migr_table_counter > 0 ? proposal->migr_table[proposal->migr_table_counter - 1].from : proposal->origin->pop; } } age1 = horizon; while (!coalesced) { //mm = mm2 = 0; //msta1 = mstart(pop1,proposal->numpop); //msto1 = mend(pop1,proposal->numpop); //msta2 = mstart(pop2,proposal->numpop); //msto2 = mend(pop2,proposal->numpop); //for (i = msta1; i < msto1; i++) //{ //mm += proposal->param0[i]; //} if(proposal->time > 1000000) return 0; mm = proposal->mig0list[pop1] * invrate; mm2 = proposal->mig0list[pop2] * invrate; if (pop1 == pop2) { denom = 2. * mm + (2. / proposal->param0[pop1] * rate); // denom_invrate = denom * invrate; proposal->time = age1 - LOG (UNIF_RANDUM ()) / denom; age1 = proposal->time; rr = UNIF_RANDUM (); r0 = ((2. / proposal->param0[pop1] * rate) / denom); r1 = r0 + mm / denom; if (rr < r0) { return 1; } else { if (rr < r1) { choice = TRUE; } else { choice = FALSE; } } } else { /*pop1 not equal pop2 */ //for (i = msta2; i < msta2; i++) //{ //mm2 += proposal->param0[i]; //} denom = mm + mm2; // denom_invrate = denom * invrate; proposal->time = age1 - LOG (UNIF_RANDUM ()) / denom; age1 = proposal->time; if (RANDUM () < (mm / denom)) { choice = TRUE; } else { choice = FALSE; } } if (choice) { if (!migrate (proposal, proposal->origin)) // &proposal->old_migr_table_counter, MIGRATION_AIR)) { return 0; /* migration limit reached */ } pop1 = proposal->migr_table_counter > 0 ? proposal->migr_table[proposal->migr_table_counter - 1].from : proposal->origin->pop; } else { if (!migrateb (proposal, proposal->realtarget)) // &proposal->old_migr_table_counter2, // MIGRATION_IN_TREE)) { return 0; /* migration limit reached */ } pop2 = proposal->migr_table_counter2 > 0 ? proposal->migr_table2[proposal->migr_table_counter2 - 1].from : proposal->realtarget->pop; } } error ("Reached the end of function without coalescing"); return -1; /*makes the compiler happy */ } #ifdef TESTING2 /// /// freeing the proposal structure, this will be replaced by a function that resets permanently /// allocated structure [reset_proposal()] void free_proposal (proposal_fmt * proposal) { // long j; long listsize = proposal->listsize; long sumtips = proposal->sumtips; long smax = proposal->world->data->maxalleles[proposal->world->locus]; // ..line_t are also reset with this memset(proposal->nodedata,0,sizeof(node *) * (2 * listsize + 2 * sumtips)); memset(proposal->mf,0, sizeof(MYREAL)*(2 * proposal->endsite)); // proposal->mt is also freed with this if (strchr (SEQUENCETYPES, proposal->datatype)) { zero_xseq(&proposal->xf,proposal->world->data->seq->endsite, proposal->world->options->rcategs); zero_xseq(&proposal->xt,proposal->world->data->seq->endsite, proposal->world->options->rcategs); } else { memset(proposal->xf.a, 0, smax); memset(proposal->xt.a, 0, smax); } memset(proposal->migr_table, 0, sizeof(migr_table_fmt) * proposal->old_migr_table_counter); memset(proposal->migr_table2, 0, sizeof(migr_table_fmt) * proposal->old_migr_table_counter2); #ifdef UEP if (proposal->world->options->uep) { memset(proposal->ueplike, 0, sizeof(MYREAL) * world->numpop * world->data->uepsites); for (j = 1; j < world->data->uepsites; ++j) (*proposal)->ueplike[j] = (*proposal)->ueplike[0] + j * world->numpop; memset(proposal->uf.s, 0, world->data->uepsites * sizeof(pair)); memset(proposal->ut.s, 0, world->data->uepsites * sizeof(pair)); memset(proposal->umf, 0, world->data->uepsites * sizeof(MYREAL)); memset(proposal->umt, 0, world->data->uepsites * sizeof(MYREAL)); } #endif // myfree(proposal); } #else /// /// freeing the proposal structure, this will be replaced by a function that resets permanently /// allocated structure [reset_proposal()] void free_proposal (proposal_fmt * proposal) { //static long count=0; myfree(proposal->aboveorigin); // myfree(proposal->bordernodes); // myfree(proposal->line_f); // myfree(proposal->line_t); // printf("%li ",count++);fflush(stdout); myfree(proposal->nodedata); // ..line_t are also freed with this myfree(proposal->mf); // proposal->mt is also freed with this if (strchr (SEQUENCETYPES, proposal->datatype)) { myfree(proposal->xf.s[0]); myfree(proposal->xt.s[0]); myfree(proposal->xf.s); myfree(proposal->xt.s); } else { myfree(proposal->xf.a); myfree(proposal->xt.a); } myfree(proposal->migr_table); myfree(proposal->migr_table2); #ifdef UEP if (proposal->world->options->uep) { myfree(proposal->ueplike[0]); myfree(proposal->ueplike); myfree(proposal->uf.s); myfree(proposal->ut.s); myfree(proposal->umf); myfree(proposal->umt); } #endif myfree(proposal); } #endif /*TESTING*/ void free_timevector (timelist_fmt * timevector) { // long i; // for (i = 0; i < timevector->allocT; i++) // { // myfree(timevector->tl[i].lineages); // } myfree(timevector->lineages); myfree(timevector->tl); myfree(timevector); } /// /// using a global timevector to save malloc/free pair calls and /// so to save time void free_timevector_new (timelist_fmt * timevector) { // memset(timevector->lineages, 0, sizeof(long) * ); // (timevector->tl); } /*----------------------------------------------------------* * rejection/acceptance of the new tree according to the likelihood * and an acceptance ratio which is higher the better the * likelihood values are (-> Metropolis) */ boolean acceptlike (world_fmt * world, proposal_fmt * proposal, long g, timelist_fmt * tyme) { //proposal changes when chain is stuck, to accept some new trees without //considering data to remove out of sticky area [hopefully] static long not_accepted = 0; #ifndef TESTINGDATE static long disturb = LONG_MAX; const long disturb_limit = proposal->world->sumtips; #endif // #ifdef UEP node *first; #endif //boolean report = TRUE; //long oldg = -1; //MYREAL zz; const long limit = MIGRATION_LIMIT * world->numpop; MYREAL rr, expo; long rm = 0L; long rmc = rmigrcount (proposal); //zz = world->heat; rm = proposal->migr_table_counter + proposal->migr_table_counter2 + world->migration_counts - rmc; if (rm > limit) { // warning might disturb parallel runs and might be not really useful at all // if (report || g != oldg) // { // warning ("migration limit (%li) exceeded: %li\n", // MIGRATION_LIMIT * world->numpop, rm); // warning // ("results may be underestimating migration rates for this chain\n"); // report = FALSE; // oldg = g; // } return FALSE; } #ifdef UEP if (world->options->uep) { first = first_uep2 (proposal->world->root->next->back, proposal->world->root->next->back, proposal->world->data->uepsites); proposal->firstuep = first_uep (first, proposal->world->root, proposal->world->data->uepsites); proposal->ueplikelihood = pseudo_ueplikelihood (world, proposal); proposal->likelihood = pseudotreelikelihood (world, proposal); } else { proposal->likelihood = pseudotreelikelihood (world, proposal); } #else // printf("DEBUG: world->likellihood[%li]=%f heat=%f\n",g,world->likelihood[g],world->heat); proposal->likelihood = pseudotreelikelihood (world, proposal); //printf("DEBUG: proposal->likellihood=%f heat=%f\n",proposal->likelihood, world->heat); #endif if (world->likelihood[g] < proposal->likelihood) { not_accepted = 0; return TRUE; } if (!world->options->heating) { expo = proposal->likelihood - world->likelihood[g]; rr = LOG (RANDUM ()); if (rr < expo) { not_accepted = 0; return TRUE; } } else { expo = (proposal->likelihood - world->likelihood[g]) * world->heat; rr = LOG (RANDUM ()); if (rr < expo) { not_accepted = 0; return TRUE; } } //printf("debug: always true\n"); //return TRUE; #ifndef TESTINGDATE if(world->cold) { if(not_accepted++ > 10000) { disturb = 0; not_accepted = 0; warning("Escaping a sticky spot by accepting the next %li genealogy changes\n (Log(P(D|G))= (previous=%f, new=%f))\n", disturb_limit,world->likelihood[g], proposal->likelihood); } if(disturb < disturb_limit) { disturb++; fprintf(stdout,"%i> %li: Log(L)=%f\n",myID, disturb, proposal->likelihood); return TRUE; } else { return FALSE; } } else { return FALSE; } #else // printf("%i> %f\t %f\t %f\n",myID, world->heat, world->likelihood[g],proposal->likelihood); return FALSE; #endif } long rmigrcount (proposal_fmt * proposal) { node *p; long count = 0; for (p = proposal->origin; p != proposal->oback; p = showtop (p->back)) { if (p->type == 'm') count++; } return count; } /// /// generates the time interval to the next event (migration or coalescence) /// and also decides what event happens. /// This funcion is modified by the rate of the mutation rate that can be /// estimated in a Bayesian context, and be fixed in a ML context /// this rate is only influencing the time and not the eventtype /// because the rate cancels out of the ratio that chooses the event. MYREAL eventtime (proposal_fmt * proposal, long pop, vtlist * tentry, char *event) { // static boolean mig_force = TRUE; MYREAL interval; MYREAL lines; MYREAL denom; MYREAL invdenom; MYREAL rate = proposal->world->options->mu_rates[proposal->world->locus]; MYREAL invrate = 1./ rate; MYREAL mm = proposal->mig0list[pop] * invrate ; //long design0 = proposal->design0list[pop]; lines = 2.0 * tentry->lineages[pop]; denom = mm + (lines * (1.0 / (proposal->param0[pop]*rate))); invdenom = 1.0 / denom; interval = (-(LOG (UNIF_RANDUM ())) * invdenom) ; //[overcautious] //if (interval < 0.0) // error("interval time is negative"); if (lines > 0.0) { if ((UNIF_RANDUM ()) < (mm * invdenom)) { *event = 'm'; return interval; } else { *event = 'c'; return interval; } } else { // printf("mcmc1.c 1653 pop = %li, lines = %li\n", pop, tentry->lineages[pop]); *event = 'm'; return interval; } } /*--------------------------------------------------------* * showsister() * find the sisternode, by going down the branch and up on * the other side again, neglecting the migration nodes. */ node * showsister (node * theNode) { node *tmp = crawlback (theNode); if (tmp->next->top) { return crawlback (tmp->next->next); } else { if (tmp->next->next->top) { return crawlback (tmp->next); } else { error ("error in treestructure, cannot find sisternode\n"); } } return NULL; } void count_migrations (node * p, long *count) { if (p->type != 't') { if (p->type == 'm') { *count += 1; count_migrations (p->next->back, count); } else { count_migrations (p->next->back, count); count_migrations (p->next->next->back, count); } } } MYREAL prob_tree (world_fmt * world, timelist_fmt * tyme) { long j, pop; MYREAL mm, cc, ss = 0; tyme->tl[0].interval = tyme->tl[0].age; for (j = 1; j < tyme->T; j++) { tyme->tl[j].interval = tyme->tl[j].age - tyme->tl[j - 1].age; } for (j = 0; j < tyme->T - 1; j++) { mm = cc = 0.0; for (pop = 0; pop < world->numpop; pop++) { mm += tyme->tl[j].lineages[pop] * world->mig0list[pop]; cc += tyme->tl[j].lineages[pop] * (tyme->tl[j].lineages[pop] - 1) / world->param0[pop]; } ss += -(tyme->tl[j].interval) * (mm + cc); if (tyme->tl[j].from == tyme->tl[j].to) ss += LOG2 - LOG (world->param0[tyme->tl[j].from]); else ss += LOG (world->param0[tyme->tl[j].from]); } return ss; }