fixed a bug with ancient samples + sweeps

Author: andrewkern

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(rn<currentFreq){
+					nodes[i]->sweepPopn=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;i<nSites;i++)printf("%d",aNode->ancSites[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;

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);

discoal_multipop.c

@@ -131,9 +131,9 @@ int main(int argc, const char * argv[]){
 
 				currentTime = sweepPhaseEventsConditionalTrajectory(&breakPoints[0], currentTime, nextTime, sweepSite, \
 					 currentFreq, &currentFreq, &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);