diff --git a/discoalFunctions.c b/discoalFunctions.c index 9870418..38b25b8 100644 --- a/discoalFunctions.c +++ b/discoalFunctions.c @@ -2241,13 +2241,20 @@ void admixPopns(int popnSrc, int popnDest1, int popnDest2, double admixProp){ } //addAncientSample -- adds ancient samples by basically flipping population id of already alloced alleles and sets there time -void addAncientSample(int lineageNumber, int popnDest, double addTime){ +void addAncientSample(int lineageNumber, int popnDest, double addTime, int stillSweeping, double currentFreq){ int i; int count = 0; + double rn; for(i=0; i < alleleNumber && count < lineageNumber; i++){ if(nodes[i]->population == (popnDest+1) * -1){ nodes[i]->population = popnDest; nodes[i]->time = addTime; + if(stillSweeping == 1){ + rn = ranf(); + if(rnsweepPopn=1; + } + } //printf("time for %d: %f\n", i, nodes[i]->time); popnSizes[popnDest]++; count++; @@ -2358,8 +2365,8 @@ rootedNode *pickNodePopnSweep(int popn,int sp){ void printNode(rootedNode *aNode){ int i; - printf("node: %p time: %f lLim: %d rLim: %d nancSites: %d popn: %d\n nsites:\n",aNode, aNode->time,aNode->lLim,\ - aNode->rLim, aNode->nancSites, aNode->population ); + printf("node: %p time: %f lLim: %d rLim: %d nancSites: %d popn: %d sweepPopn: %d\n",aNode, aNode->time,aNode->lLim,\ + aNode->rLim, aNode->nancSites, aNode->population, aNode->sweepPopn); for(i=0;iancSites[i]); printf("\n"); } @@ -2409,6 +2416,13 @@ void printAllNodes(){ } } +void printAllActiveNodes(){ + int i; + for(i = 0 ; i < alleleNumber; i++){ + printNode(nodes[i]); + } +} + /********** event stuff */ int compare_events(const void *a,const void *b){ struct event *pa = (struct event *) a; diff --git a/discoalFunctions.h b/discoalFunctions.h index 57ea1ac..930dd28 100644 --- a/discoalFunctions.h +++ b/discoalFunctions.h @@ -30,7 +30,7 @@ void errorCheckMutations(); void mergePopns(int popnSrc, int popnDest); void admixPopns(int popnSrc, int popnDest1, int popnDest2, double admixProp); -void addAncientSample(int lineageNumber, int popnDest, double addTime); +void addAncientSample(int lineageNumber, int popnDest, double addTime, int stillSweeping, double currentFreq); void recurrentMutAtTime(double cTime,int srcPopn, int sp); @@ -83,6 +83,7 @@ int isCoalNode(rootedNode *aNode); void newickRecurse(rootedNode *aNode, float site,float tempTime); void printTreeAtSite(float site); void printAllNodes(); +void printAllActiveNodes(); unsigned int devrand(void); int compare_doubles(const void *a,const void *b); diff --git a/discoal_multipop.c b/discoal_multipop.c index 66187b4..0643b2b 100644 --- a/discoal_multipop.c +++ b/discoal_multipop.c @@ -131,9 +131,9 @@ int main(int argc, const char * argv[]){ currentTime = sweepPhaseEventsConditionalTrajectory(&breakPoints[0], currentTime, nextTime, sweepSite, \ currentFreq, ¤tFreq, &activeSweepFlag, alpha, currentSize, sweepMode, f0, uA); - // printf("currentFreqAfter: %f alleleNumber:%d\n",currentFreq,alleleNumber); - // printf("pn0:%d pn1:%d alleleNumber: %d sp1: %d sp2: %d \n", popnSizes[0],popnSizes[1], alleleNumber,sweepPopnSizes[1], - // sweepPopnSizes[0]); + //printf("currentFreqAfter: %f alleleNumber:%d currentTime:%f\n",currentFreq,alleleNumber,currentTime); + //printf("pn0:%d pn1:%d alleleNumber: %d sp1: %d sp2: %d \n", popnSizes[0],popnSizes[1], alleleNumber,sweepPopnSizes[1], + // sweepPopnSizes[0]); if (currentTime < nextTime) currentTime = neutralPhaseGeneralPopNumber(&breakPoints[0], currentTime, nextTime, currentSize); @@ -194,7 +194,9 @@ int main(int argc, const char * argv[]){ break; case 'A': currentTime = events[j].time; - addAncientSample(events[j].lineageNumber, events[j].popID, events[j].time); + //printAllActiveNodes(); + addAncientSample(events[j].lineageNumber, events[j].popID, events[j].time, activeSweepFlag, currentFreq); + //printAllActiveNodes(); if(activeSweepFlag == 0){ if(recurSweepMode ==0){ currentTime = neutralPhaseGeneralPopNumber(&breakPoints[0], currentTime, nextTime, currentSize);