##########################function definition ###########################3


#function to estimate the heterozygosity of a sample

het.funct<-function(sample){

	nall<-ncol(sample)
	nloc<-nall/2
	IndHetLoci<-1:nrow(sample)

	for(ind in 1:nrow(sample))
	IndHetLoci[ind]<-length(which( sample[ind,seq(1,nall-1,2)] != sample[ind,seq(2,nall,2)]   ))/nloc

	IndHetLoci
}

########################## end of function definition ###########################3


#This reads the appropriate sample file into R (Just delete the #, or add a new file name)
#Reading the file into R, defining it as "parents"
parents<-read.delim("Carp021612RNoMpy9.txt",header=FALSE)
#parents<-read.delim("Gol021612RNoMpy9.txt",header=FALSE)
#parents<-read.delim("Moh021612RNoMpy9.txt",header=FALSE)

nrep<-10000

PTested<-seq(0,1,0.01)

Cull0<-matrix(nrow=length(PTested),ncol=4)
Cull1<-matrix(nrow=length(PTested),ncol=4)
Cull2<-matrix(nrow=length(PTested),ncol=4)
Cull3<-matrix(nrow=length(PTested),ncol=4)
Cull4<-matrix(nrow=length(PTested),ncol=4)


CullR<-1

### SELFING, NO MORTALITY ###
for(Prop in PTested){
	Selfed<-nrep*Prop

	#Defining the number of alleles as the number of columns in the file
	nall<-ncol(parents)
	nloc<-nall/2
	#Defining the number of individuals as the number of rows in the "parents" file
	nind<-nrow(parents)

	#Random index, samples from the "parents" file to give a random sampling of parents which will be
	#used to generate our multilocus genotypes
	Rindex1<-sample(1:nind, Selfed+1, replace=T)
	Rindex2<-sample(1:nind, (2*(nrep-Selfed)), replace=T)


	#Creates a matrix of the random parental genotypes, based on the Rindex above.
	Rparents1<-parents[Rindex1,]
	Rparents2<-parents[Rindex2,]

	#Creates a matrix to hold our single offspring multilocus genotypes
	MLGone<-matrix(nrow=nrep,ncol=nall)


	#This for loop randomly samples the parental genotypes in Rparents to create 1 offspring genotype per parent.
	#It then places the offspring genotypes in the MLGpair matrix
	if(Selfed>0){
		for(ind in 1:(Selfed)){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGone[ind,start:stop]<-sort(sample(as.numeric(Rparents1[ind,start:stop]),2,replace=T))
			}
		}
	}
		
	Selfed1<-Selfed+1
	for(ind in seq(1,(2*nrep),2)){
		if(Selfed1<=nrep){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGone[Selfed1,start]<-sample(as.numeric(Rparents2[ind,start:stop]),1)
				MLGone[Selfed1,stop]<-sample(as.numeric(Rparents2[(ind+1),start:stop]),1)
				MLGone[Selfed1,start:stop]<-sort(as.numeric(MLGone[Selfed1,start:stop]))
			}
			Selfed1<-Selfed1+1
		}	
	}

	#Important Histogram Definitions
	#Defines H0SelfedHet as the result of our het.funct (above) as implemented on the MLGone matrix (simulated data)
	H0SelfedHet<-het.funct(MLGone)
	#Defines Observed.Het as the result of our het.funct (above) as implemented on the original data
	Observed.Het<-het.funct(parents)
	
	nH0<-length(H0SelfedHet)
	nObs<-length(Observed.Het)
	
		
	Proportions<-matrix(nrow=11,ncol=7)
	
	class<-round(0,1)
	for(r in seq(1,11,1)){
		# Number per class (Simulated)
		Proportions[r,1]<-sum(H0SelfedHet==(class))
		# Number per class (Observed)
		Proportions[r,2]<-sum(Observed.Het==(class))
		# Proportion per class (Simulated)
		Proportions[r,3]<-Proportions[r,1]/nH0
		# Proportion per class (Observed) 
		Proportions[r,4]<-Proportions[r,2]/nObs
		# Proportion difference
		Proportions[r,5]<-(Proportions[r,3]-Proportions[r,4])
		# Corrected Proportion Difference (Modified Below)
		Proportions[r,6]<-Proportions[r,5]
		class<-round(class+0.1,1)
	}
	
	# For loop to correct proportion differences (eliminates negatives and differences for high No. Het Loci)
	for(r in seq(1,11,1)){
		if (Proportions[r,6]<0){
			Proportions[r,6]<-0
		} 
		if (r > 5){
			if (Proportions[r-1,6]==0){
				if (Proportions[r,6]>0){
					Proportions[r,6]<-0
				}
			}
		}
	}
	
	# Calculates Proportion Surviving, Places into Column 7	
	Proportions[1:11,7]<-1 ###Altered to eliminate mortality###
	#(1-(Proportions[1:11,6]/Proportions[1:11,3]))
	
	# Eliminates NaN from matrix (NaN occur when dividing by zero)
	mask<-apply(Proportions, 2, is.nan)
	Proportions[mask]<-1
	
	# Places sum of Proportion Culled into Cull matrix
	Cull0[CullR,1]<-sum(Proportions[1:11,6])
	
		
	####################   Mortality   ####################

	Mort<-Proportions[1:11,1]*Proportions[1:11,7]

	MortHet<-round(rep(seq(0,1,0.1), times=Mort),1)
	nMort<-length(MortHet)

	MortProp<-matrix(nrow=11,ncol=7)
	
	class<-round(0,1)
	for(r in seq(1,11,1)){
	
		# Number per class (Simulated, Post-Mortality)
		MortProp[r,1]<-sum(MortHet==(class))
		# Number per class (Observed)
		MortProp[r,2]<-sum(Observed.Het==(class))
		# Proportion per class (Simulated, Post-Mortality)
		MortProp[r,3]<-MortProp[r,1]/nMort
		# Proportion per class (Observed) 
		MortProp[r,4]<-MortProp[r,2]/nObs
		# Absolute Proportion difference
		MortProp[r,5]<-sqrt((MortProp[r,4]-MortProp[r,3])^2)
		# Number Simulated per class
		MortProp[r,6]<-MortProp[r,3]*nObs
		if (MortProp[r,6] < 0.01){
			MortProp[r,6]<-0.01
		}
		# Chi2 Per Class
		MortProp[r,7]<-((MortProp[r,2]-MortProp[r,6])^2)/(MortProp[r,6])
		
						
		class<-round(class+0.1,1)
	}
	
	# Eliminates NaN from matrix (NaN occur when dividing by zero)
	mask<-apply(MortProp, 2, is.nan)
	MortProp[mask]<-0
	
	
	Cull0[CullR,2]<-sum(MortProp[1:11,5])
	
	Cull0[CullR,3]<-sum(MortProp[1:11,7])
	
	Cull0[CullR,4]<-(1-pchisq(Cull0[CullR,3],nrow(MortProp)-1))
	
	CullR<-CullR+1
	
	print(Prop)
}

### SELFING, WITH MORTALITY ###
CullR<-1

for(Prop in PTested){
	Selfed<-nrep*Prop

	#Defining the number of alleles as the number of columns in the file
	nall<-ncol(parents)
	nloc<-nall/2
	#Defining the number of individuals as the number of rows in the "parents" file
	nind<-nrow(parents)

	#Random index, samples from the "parents" file to give a random sampling of parents which will be
	#used to generate our multilocus genotypes
	Rindex1<-sample(1:nind, Selfed+1, replace=T)
	Rindex2<-sample(1:nind, (2*(nrep-Selfed)), replace=T)


	#Creates a matrix of the random parental genotypes, based on the Rindex above.
	Rparents1<-parents[Rindex1,]
	Rparents2<-parents[Rindex2,]

	#Creates a matrix to hold our single offspring multilocus genotypes
	MLGone<-matrix(nrow=nrep,ncol=nall)


	#This for loop randomly samples the parental genotypes in Rparents to create 1 offspring genotype per parent.
	#It then places the offspring genotypes in the MLGpair matrix
	if(Selfed>0){
		for(ind in 1:(Selfed)){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGone[ind,start:stop]<-sort(sample(as.numeric(Rparents1[ind,start:stop]),2,replace=T))
			}
		}
	}
		
	Selfed1<-Selfed+1
	for(ind in seq(1,(2*nrep),2)){
		if(Selfed1<=nrep){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGone[Selfed1,start]<-sample(as.numeric(Rparents2[ind,start:stop]),1)
				MLGone[Selfed1,stop]<-sample(as.numeric(Rparents2[(ind+1),start:stop]),1)
				MLGone[Selfed1,start:stop]<-sort(as.numeric(MLGone[Selfed1,start:stop]))
			}
			Selfed1<-Selfed1+1
		}	
	}

	#Important Histogram Definitions
	#Defines H0SelfedHet as the result of our het.funct (above) as implemented on the MLGone matrix (simulated data)
	H0SelfedHet<-het.funct(MLGone)
	#Defines Observed.Het as the result of our het.funct (above) as implemented on the original data
	Observed.Het<-het.funct(parents)
	
	nH0<-length(H0SelfedHet)
	nObs<-length(Observed.Het)
	
		
	Proportions<-matrix(nrow=11,ncol=7)
	
	class<-round(0,1)
	for(r in seq(1,11,1)){
		# Number per class (Simulated)
		Proportions[r,1]<-sum(H0SelfedHet==(class))
		# Number per class (Observed)
		Proportions[r,2]<-sum(Observed.Het==(class))
		# Proportion per class (Simulated)
		Proportions[r,3]<-Proportions[r,1]/nH0
		# Proportion per class (Observed) 
		Proportions[r,4]<-Proportions[r,2]/nObs
		# Proportion difference
		Proportions[r,5]<-(Proportions[r,3]-Proportions[r,4])
		# Corrected Proportion Difference (Modified Below)
		Proportions[r,6]<-Proportions[r,5]
		class<-round(class+0.1,1)
	}
	
	# For loop to correct proportion differences (eliminates negatives and differences for high No. Het Loci)
	for(r in seq(1,11,1)){
		if (Proportions[r,6]<0){
			Proportions[r,6]<-0
		} 
		if (r > 5){
			if (Proportions[r-1,6]==0){
				if (Proportions[r,6]>0){
					Proportions[r,6]<-0
				}
			}
		}
	}
	
	# Calculates Proportion Surviving, Places into Column 7	
	Proportions[1:11,7]<-(1-(Proportions[1:11,6]/Proportions[1:11,3]))
	
	# Eliminates NaN from matrix (NaN occur when dividing by zero)
	mask<-apply(Proportions, 2, is.nan)
	Proportions[mask]<-1
	
	# Places sum of Proportion Culled into Cull matrix
	Cull1[CullR,1]<-sum(Proportions[1:11,6])
	
		
	####################   Mortality   ####################

	Mort<-Proportions[1:11,1]*Proportions[1:11,7]

	MortHet<-round(rep(seq(0,1,0.1), times=Mort),1)
	nMort<-length(MortHet)

	MortProp<-matrix(nrow=11,ncol=7)
	
	class<-round(0,1)
	for(r in seq(1,11,1)){
	
		# Number per class (Simulated, Post-Mortality)
		MortProp[r,1]<-sum(MortHet==(class))
		# Number per class (Observed)
		MortProp[r,2]<-sum(Observed.Het==(class))
		# Proportion per class (Simulated, Post-Mortality)
		MortProp[r,3]<-MortProp[r,1]/nMort
		# Proportion per class (Observed) 
		MortProp[r,4]<-MortProp[r,2]/nObs
		# Absolute Proportion difference
		MortProp[r,5]<-sqrt((MortProp[r,4]-MortProp[r,3])^2)
		# Number Simulated per class
		MortProp[r,6]<-MortProp[r,3]*nObs
		if (MortProp[r,6] < 0.01){
			MortProp[r,6]<-0.01
		}
		# Chi2 Per Class
		MortProp[r,7]<-((MortProp[r,2]-MortProp[r,6])^2)/(MortProp[r,6])
		
						
		class<-round(class+0.1,1)
	}
	
	# Eliminates NaN from matrix (NaN occur when dividing by zero)
	mask<-apply(MortProp, 2, is.nan)
	MortProp[mask]<-0
	
	
	Cull1[CullR,2]<-sum(MortProp[1:11,5])
	
	Cull1[CullR,3]<-sum(MortProp[1:11,7])
	
	Cull1[CullR,4]<-(1-pchisq(Cull1[CullR,3],nrow(MortProp)-1))
	
	CullR<-CullR+1
	
	print(Prop)
}

### FULL SIBLINGS ###
CullR<-1

for(PropSibs in PTested){
	#Defining the number of alleles as the number of columns in the file
	nall<-ncol(parents)
	nloc<-nall/2
	Sibs<-nrep*PropSibs

	#Defining the number of individuals as the number of rows in the "parents" file
	nind<-nrow(parents)

	#Random index, samples from the "parents" file to give a random sampling of parents which will be
	#used to generate our multilocus genotypes
	Rindex1<-sample(1:nind, Sibs, replace=T)
	Rindex2<-sample(1:nind, (2*(nrep-Sibs)), replace=T)


	#Creates a matrix of the random parental genotypes, based on the Rindex above.
	Rparents1<-parents[Rindex1,]
	Rparents2<-parents[Rindex2,]

	#Creates a matrix to hold our single offspring multilocus genotypes
	MLGSibs<-matrix(nrow=nrep,ncol=nall)


	#This for loop randomly samples the parental genotypes in Rparents1 to create 2 Full-sib genotypes.
	#It then places the offspring genotypes in the MLGSibs matrix
	if(PropSibs>0){
		for(ind in seq(1,Sibs,2)){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGSibs[ind:(ind+1),start]<-sample(as.numeric(Rparents1[ind,start:stop]),2, replace=T)
				MLGSibs[ind:(ind+1),stop]<-sample(as.numeric(Rparents1[(ind+1),start:stop]),2, replace=T)
				MLGSibs[ind:(ind+1),start:stop]<-sort(as.numeric(MLGSibs[ind:(ind+1),start:stop]))
			}
		}
	}

	#This for loop randomly samples parental genotypes in Rparents2 to create outcrossed offspring
	Sibs1<-Sibs+1
	for(ind in seq(1,(2*nrep),2)){
		if(Sibs1<(nrep+1)){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGSibs[Sibs1,start]<-sample(as.numeric(Rparents2[ind,start:stop]),1)
				MLGSibs[Sibs1,stop]<-sample(as.numeric(Rparents2[(ind+1),start:stop]),1)
				MLGSibs[Sibs1,start:stop]<-sort(as.numeric(MLGSibs[Sibs1,start:stop]))
			}
			Sibs1<-Sibs1+1
		}	
	}

	#Defines H0SelfedHet as the result of our het.funct (above) as implemented on the MLGone matrix (simulated data)
	H0SelfedHet<-het.funct(MLGSibs)
	#Defines Observed.Het as the result of our het.funct (above) as implemented on the original data
	Observed.Het<-het.funct(parents)
	
	nH0<-length(H0SelfedHet)
	nObs<-length(Observed.Het)
	
	Proportions<-matrix(nrow=11,ncol=7)
	
	class<-round(0,1)
	for(r in seq(1,11,1)){
		# Number per class (Simulated)
		Proportions[r,1]<-sum(H0SelfedHet==(class))
		# Number per class (Observed)
		Proportions[r,2]<-sum(Observed.Het==(class))
		# Proportion per class (Simulated)
		Proportions[r,3]<-Proportions[r,1]/nH0
		# Proportion per class (Observed) 
		Proportions[r,4]<-Proportions[r,2]/nObs
		# Proportion difference
		Proportions[r,5]<-(Proportions[r,3]-Proportions[r,4])
		# Corrected Proportion Difference (Modified Below)
		Proportions[r,6]<-Proportions[r,5]
		class<-round(class+0.1,1)
	}
	
	# For loop to correct proportion differences (eliminates negatives and differences for high No. Het Loci)
	for(r in seq(1,11,1)){
		if (Proportions[r,6]<0){
			Proportions[r,6]<-0
		} 
		if (r > 5){
			if (Proportions[r-1,6]==0){
				if (Proportions[r,6]>0){
					Proportions[r,6]<-0
				}
			}
		}
	}
	
	# Calculates Proportion Surviving, Places into Column 7	
	Proportions[1:11,7]<-(1-(Proportions[1:11,6]/Proportions[1:11,3]))
	
	# Eliminates NaN from matrix (NaN occur when dividing by zero)
	mask<-apply(Proportions, 2, is.nan)
	Proportions[mask]<-1
	
	# Places sum of Proportion Culled into Cull matrix
	Cull2[CullR,1]<-sum(Proportions[1:11,6])
	
		
	####################   Mortality   ####################

	Mort<-Proportions[1:11,1]*Proportions[1:11,7]

	MortHet<-round(rep(seq(0,1,0.1), times=Mort),1)
	nMort<-length(MortHet)

	MortProp<-matrix(nrow=11,ncol=7)
	
	class<-round(0,1)
	for(r in seq(1,11,1)){
	
		# Number per class (Simulated, Post-Mortality)
		MortProp[r,1]<-sum(MortHet==(class))
		# Number per class (Observed)
		MortProp[r,2]<-sum(Observed.Het==(class))
		# Proportion per class (Simulated, Post-Mortality)
		MortProp[r,3]<-MortProp[r,1]/nMort
		# Proportion per class (Observed) 
		MortProp[r,4]<-MortProp[r,2]/nObs
		# Absolute Proportion difference
		MortProp[r,5]<-sqrt((MortProp[r,4]-MortProp[r,3])^2)
		# Number Simulated per class
		MortProp[r,6]<-MortProp[r,3]*nObs
		if (MortProp[r,6] < 0.01){
			MortProp[r,6]<-0.01
		}
		# Chi2 Per Class
		MortProp[r,7]<-((MortProp[r,2]-MortProp[r,6])^2)/(MortProp[r,6])
		
						
		class<-round(class+0.1,1)
	}
	
	# Eliminates NaN from matrix (NaN occur when dividing by zero)
	mask<-apply(MortProp, 2, is.nan)
	MortProp[mask]<-0
	
	
	Cull2[CullR,2]<-sum(MortProp[1:11,5])
	
	Cull2[CullR,3]<-sum(MortProp[1:11,7])
	
	Cull2[CullR,4]<-(1-pchisq(Cull2[CullR,3],nrow(MortProp)-1))
	
	CullR<-CullR+1
	
	print(PropSibs)
}

### HALF SIBS ###
CullR<-1


for(PropSibs in PTested){
	#Defining the number of alleles as the number of columns in the file
	nall<-ncol(parents)
	nloc<-nall/2
	Sibs<-nrep*PropSibs

	#Defining the number of individuals as the number of rows in the "parents" file
	nind<-nrow(parents)

	#Random index, samples from the "parents" file to give a random sampling of parents which will be
	#used to generate our multilocus genotypes
	Rindex1<-sample(1:nind, ((3*Sibs)/2), replace=T)
	Rindex2<-sample(1:nind, (2*(nrep-Sibs)), replace=T)


	#Creates a matrix of the random parental genotypes, based on the Rindex above.
	Rparents1<-parents[Rindex1,]
	Rparents2<-parents[Rindex2,]

	#Creates a matrix to hold our single offspring multilocus genotypes
	MLGSibs<-matrix(nrow=nrep,ncol=nall)


	#This for loop randomly samples the parental genotypes in Rparents1 to create 2 Half-sib genotypes.
	#It then places the offspring genotypes in the MLGSibs matrix
	par<-1
	if(PropSibs>0){
		for(ind in seq(1,Sibs,2)){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGSibs[ind:(ind+1),start]<-sample(as.numeric(Rparents1[par,start:stop]),2, replace=T)
				MLGSibs[ind,stop]<-sample(as.numeric(Rparents1[(par+1),start:stop]),1, replace=T)
				MLGSibs[(ind+1),stop]<-sample(as.numeric(Rparents1[(par+2),start:stop]),1, replace=T)
				MLGSibs[ind:(ind+1),start:stop]<-sort(as.numeric(MLGSibs[ind:(ind+1),start:stop]))
			}
			par<-par+3
		}
	}

	#This for loop randomly samples parental genotypes in Rparents2 to create outcrossed offspring
	Sibs1<-Sibs+1
	for(ind in seq(1,(2*nrep),2)){
		if(Sibs1<(nrep+1)){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGSibs[Sibs1,start]<-sample(as.numeric(Rparents2[ind,start:stop]),1)
				MLGSibs[Sibs1,stop]<-sample(as.numeric(Rparents2[(ind+1),start:stop]),1)
				MLGSibs[Sibs1,start:stop]<-sort(as.numeric(MLGSibs[Sibs1,start:stop]))
			}
			Sibs1<-Sibs1+1
		}	
	}

	#Important Histogram Definitions
	#Defines H0SelfedHet as the result of our het.funct (above) as implemented on the MLGone matrix (simulated data)
	H0SelfedHet<-het.funct(MLGSibs)
	
	#Defines Observed.Het as the result of our het.funct (above) as implemented on the original data
	Observed.Het<-het.funct(parents)
	
	nH0<-length(H0SelfedHet)
	nObs<-length(Observed.Het)

	
	Proportions<-matrix(nrow=11,ncol=7)
	
	class<-round(0,1)
	for(r in seq(1,11,1)){
		# Number per class (Simulated)
		Proportions[r,1]<-sum(H0SelfedHet==(class))
		# Number per class (Observed)
		Proportions[r,2]<-sum(Observed.Het==(class))
		# Proportion per class (Simulated)
		Proportions[r,3]<-Proportions[r,1]/nH0
		# Proportion per class (Observed) 
		Proportions[r,4]<-Proportions[r,2]/nObs
		# Proportion difference
		Proportions[r,5]<-(Proportions[r,3]-Proportions[r,4])
		# Corrected Proportion Difference (Modified Below)
		Proportions[r,6]<-Proportions[r,5]
		class<-round(class+0.1,1)
	}
	
	# For loop to correct proportion differences (eliminates negatives and differences for high No. Het Loci)
	for(r in seq(1,11,1)){
		if (Proportions[r,6]<0){
			Proportions[r,6]<-0
		} 
		if (r > 5){
			if (Proportions[r-1,6]==0){
				if (Proportions[r,6]>0){
					Proportions[r,6]<-0
				}
			}
		}
	}
	
	# Calculates Proportion Surviving, Places into Column 7	
	Proportions[1:11,7]<-(1-(Proportions[1:11,6]/Proportions[1:11,3]))
	
	# Eliminates NaN from matrix (NaN occur when dividing by zero)
	mask<-apply(Proportions, 2, is.nan)
	Proportions[mask]<-1
	
	# Places sum of Proportion Culled into Cull matrix
	Cull3[CullR,1]<-sum(Proportions[1:11,6])
	
		
	####################   Mortality   ####################

	Mort<-Proportions[1:11,1]*Proportions[1:11,7]

	MortHet<-round(rep(seq(0,1,0.1), times=Mort),1)
	nMort<-length(MortHet)

	MortProp<-matrix(nrow=11,ncol=7)
	
	class<-round(0,1)
	for(r in seq(1,11,1)){
	
		# Number per class (Simulated, Post-Mortality)
		MortProp[r,1]<-sum(MortHet==(class))
		# Number per class (Observed)
		MortProp[r,2]<-sum(Observed.Het==(class))
		# Proportion per class (Simulated, Post-Mortality)
		MortProp[r,3]<-MortProp[r,1]/nMort
		# Proportion per class (Observed) 
		MortProp[r,4]<-MortProp[r,2]/nObs
		# Absolute Proportion difference
		MortProp[r,5]<-sqrt((MortProp[r,4]-MortProp[r,3])^2)
		# Number Simulated per class
		MortProp[r,6]<-MortProp[r,3]*nObs
		if (MortProp[r,6] < 0.01){
			MortProp[r,6]<-0.01
		}
		# Chi2 Per Class
		MortProp[r,7]<-((MortProp[r,2]-MortProp[r,6])^2)/(MortProp[r,6])
		
						
		class<-round(class+0.1,1)
	}
	
	# Eliminates NaN from matrix (NaN occur when dividing by zero)
	mask<-apply(MortProp, 2, is.nan)
	MortProp[mask]<-0
	
	
	Cull3[CullR,2]<-sum(MortProp[1:11,5])
	
	Cull3[CullR,3]<-sum(MortProp[1:11,7])
	
	Cull3[CullR,4]<-(1-pchisq(Cull3[CullR,3],nrow(MortProp)-1))
	
	CullR<-CullR+1
	
	print(PropSibs)
}

### FIRST COUSINS ###
CullR<-1


for(PropCous in PTested){
	#Defining the number of alleles as the number of columns in the file
	nall<-ncol(parents)
	nloc<-nall/2
	Cous<-nrep*PropCous

	#Defining the number of individuals as the number of rows in the "parents" file
	nind<-nrow(parents)

	#Random indexes, generate random number list from 1:(user defined), used to create random lists of parents, below
	Rindex1<-sample(1:nind, Cous, replace=T)
	Rindex2<-sample(1:nind, (2*Cous), replace=T)
	Rindex3<-sample(1:nind, 2*(nrep-Cous),replace=T)

	#Creates matrices of the random parental genotypes, based on the Rindexes above.
	Rparents1<-parents[Rindex1,]
	Rparents2<-parents[Rindex2,]
	Rparents3<-parents[Rindex3,]

	#Creates matrices to hold our offspring multilocus genotypes
	MLGSibs<-matrix(nrow=(nrep*PropCous),ncol=nall)
	MLGOut<-matrix(nrow=(nrep*PropCous),ncol=nall)
	MLGCous<-matrix(nrow=nrep,ncol=nall)


	#This for loop randomly samples the parental genotypes in Rparents1 to create 2 Full-sib genotypes.
	#It then places the offspring genotypes in the MLGSibs matrix
	#start<-1
	#ind<-1
	Out1<-1
	if(PropCous>0){
		for(ind in seq(1,Cous,2)){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGSibs[ind:(ind+1),start]<-sample(as.numeric(Rparents1[ind,start:stop]),2, replace=T)
				MLGSibs[ind:(ind+1),stop]<-sample(as.numeric(Rparents1[(ind+1),start:stop]),2, replace=T)
				MLGSibs[ind:(ind+1),start:stop]<-sort(as.numeric(MLGSibs[ind:(ind+1),start:stop]))
			}
		}
		for(ind in seq(1,2*(Cous),2)){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGOut[Out1,start]<-sample(as.numeric(Rparents2[ind,start:stop]),1)
				MLGOut[Out1,stop]<-sample(as.numeric(Rparents2[(ind+1),start:stop]),1)
				MLGOut[Out1,start:stop]<-sort(as.numeric(MLGOut[Out1,start:stop]))
			}
			Out1<-Out1+1
		}
		for(ind in 1:Cous){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGCous[ind,start]<-sample(as.numeric(MLGSibs[ind,start:stop]),1)
				MLGCous[ind,stop]<-sample(as.numeric(MLGOut[ind,start:stop]),1)
				MLGCous[ind,start:stop]<-sort(as.numeric(MLGCous[ind,start:stop]))
			}
		}
	}
	
	#This for loop randomly samples parental genotypes in Rparents3 to create outcrossed offspring
	Cous1<-Cous+1
	for(ind in seq(1,(2*nrep),2)){
		if(Cous1<(nrep+1)){
			for(start in seq(1,(nall-1),2)){
				stop<-start+1
				MLGCous[Cous1,start]<-sample(as.numeric(Rparents3[ind,start:stop]),1)
				MLGCous[Cous1,stop]<-sample(as.numeric(Rparents3[(ind+1),start:stop]),1)
				MLGCous[Cous1,start:stop]<-sort(as.numeric(MLGCous[Cous1,start:stop]))
			}
			Cous1<-Cous1+1
		}	
	}

	

	#Important Histogram Definitions
	#Defines H0SelfedHet as the result of our het.funct (above) as implemented on the MLGone matrix (simulated data)
	H0SelfedHet<-het.funct(MLGCous)
	
	#Defines Observed.Het as the result of our het.funct (above) as implemented on the original data
	Observed.Het<-het.funct(parents)
	
	nH0<-length(H0SelfedHet)
	nObs<-length(Observed.Het)


	Proportions<-matrix(nrow=11,ncol=7)
	
	class<-round(0,1)
	for(r in seq(1,11,1)){
		# Number per class (Simulated)
		Proportions[r,1]<-sum(H0SelfedHet==(class))
		# Number per class (Observed)
		Proportions[r,2]<-sum(Observed.Het==(class))
		# Proportion per class (Simulated)
		Proportions[r,3]<-Proportions[r,1]/nH0
		# Proportion per class (Observed) 
		Proportions[r,4]<-Proportions[r,2]/nObs
		# Proportion difference
		Proportions[r,5]<-(Proportions[r,3]-Proportions[r,4])
		# Corrected Proportion Difference (Modified Below)
		Proportions[r,6]<-Proportions[r,5]
		class<-round(class+0.1,1)
	}
	
	# For loop to correct proportion differences (eliminates negatives and differences for high No. Het Loci)
	for(r in seq(1,11,1)){
		if (Proportions[r,6]<0){
			Proportions[r,6]<-0
		} 
		if (r > 5){
			if (Proportions[r-1,6]==0){
				if (Proportions[r,6]>0){
					Proportions[r,6]<-0
				}
			}
		}
	}
	
	# Calculates Proportion Surviving, Places into Column 7	
	Proportions[1:11,7]<-(1-(Proportions[1:11,6]/Proportions[1:11,3]))
	
	# Eliminates NaN from matrix (NaN occur when dividing by zero)
	mask<-apply(Proportions, 2, is.nan)
	Proportions[mask]<-1
	
	# Places sum of Proportion Culled into Cull matrix
	Cull4[CullR,1]<-sum(Proportions[1:11,6])
	
		
	####################   Mortality   ####################

	Mort<-Proportions[1:11,1]*Proportions[1:11,7]

	MortHet<-round(rep(seq(0,1,0.1), times=Mort),1)
	nMort<-length(MortHet)

	MortProp<-matrix(nrow=11,ncol=7)
	
	class<-round(0,1)
	for(r in seq(1,11,1)){
	
		# Number per class (Simulated, Post-Mortality)
		MortProp[r,1]<-sum(MortHet==(class))
		# Number per class (Observed)
		MortProp[r,2]<-sum(Observed.Het==(class))
		# Proportion per class (Simulated, Post-Mortality)
		MortProp[r,3]<-MortProp[r,1]/nMort
		# Proportion per class (Observed) 
		MortProp[r,4]<-MortProp[r,2]/nObs
		# Absolute Proportion difference
		MortProp[r,5]<-sqrt((MortProp[r,4]-MortProp[r,3])^2)
		# Number Simulated per class
		MortProp[r,6]<-MortProp[r,3]*nObs
		if (MortProp[r,6] < 0.01){
			MortProp[r,6]<-0.01
		}
		# Chi2 Per Class
		MortProp[r,7]<-((MortProp[r,2]-MortProp[r,6])^2)/(MortProp[r,6])
		
						
		class<-round(class+0.1,1)
	}
	
	# Eliminates NaN from matrix (NaN occur when dividing by zero)
	mask<-apply(MortProp, 2, is.nan)
	MortProp[mask]<-0
	
	
	Cull4[CullR,2]<-sum(MortProp[1:11,5])
	
	Cull4[CullR,3]<-sum(MortProp[1:11,7])
	
	Cull4[CullR,4]<-(1-pchisq(Cull4[CullR,3],nrow(MortProp)-1))
	
	CullR<-CullR+1
	
	print(PropCous)
}

#Adding a column of proportions tested
Cull0<-cbind(Cull0[,c(1:4)],(PTested*100))
Cull1<-cbind(Cull1[,c(1:4)],(PTested*100))
Cull2<-cbind(Cull2[,c(1:4)],(PTested*100))
Cull3<-cbind(Cull3[,c(1:4)],(PTested*100))
Cull4<-cbind(Cull4[,c(1:4)],(PTested*100))

#Significance range of Chi Square values
which(Cull1[,4]>=0.05)

### SAVING OUTPUT TO FILES ###

write.table(Cull0, file="CarpCull0.txt", sep="\t")
write.table(Cull1, file="CarpCull1.txt", sep="\t")
write.table(Cull2, file="CarpCull2.txt", sep="\t")
write.table(Cull3, file="CarpCull3.txt", sep="\t")
write.table(Cull4, file="CarpCull4.txt", sep="\t")

write.table(Cull0, file="GolCull0.txt", sep="\t")
write.table(Cull1, file="GolCull1.txt", sep="\t")
write.table(Cull2, file="GolCull2.txt", sep="\t")
write.table(Cull3, file="GolCull3.txt", sep="\t")
write.table(Cull4, file="GolCull4.txt", sep="\t")

write.table(Cull0, file="MohCull0.txt", sep="\t")
write.table(Cull1, file="MohCull1.txt", sep="\t")
write.table(Cull2, file="MohCull2.txt", sep="\t")
write.table(Cull3, file="MohCull3.txt", sep="\t")
write.table(Cull4, file="MohCull4.txt", sep="\t")

### READING INTO R ###
Carp0<-read.delim("CarpCull0.txt")
Carp1<-read.delim("CarpCull1.txt")
Carp2<-read.delim("CarpCull2.txt")
Carp3<-read.delim("CarpCull3.txt")
Carp4<-read.delim("CarpCull4.txt")

Gol0<-read.delim("GolCull0.txt")
Gol1<-read.delim("GolCull1.txt")
Gol2<-read.delim("GolCull2.txt")
Gol3<-read.delim("GolCull3.txt")
Gol4<-read.delim("GolCull4.txt")

Moh0<-read.delim("MohCull0.txt")
Moh1<-read.delim("MohCull1.txt")
Moh2<-read.delim("MohCull2.txt")
Moh3<-read.delim("MohCull3.txt")
Moh4<-read.delim("MohCull4.txt")


### PLOTTING FOR PAPER ###

########################   SINGLE PLOTS   ######################################
#Carpinteria
windows.options(width=7, height=10)
par(mfrow=c(3,1), bty="n", cex.axis=1.5, cex.lab=1.5)

par(mar=c(1,7,0,0), bty="n")
plot(predict(loess(Carp1[,2]~Carp1[,5],span=0.5)), col="black", type="l", lwd=3, lty=1, xlab="",
	ylab="Absolute Difference", ylim=c(0,1.0), xlim=c(0,100), axes=F)
lines(which.min(predict(loess(Carp1[,2]~Carp1[,5],span=0.5))), min(predict(loess(Carp1[,2]~Carp1[,5],span=0.5))),
	type="p", col="black", pch=5, cex=1.5)
lines(predict(loess(Carp0[,2]~Carp0[,5],span=0.5)), col="gray50", lwd=3, lty=1)
lines(predict(loess(Carp2[,2]~Carp2[,5],span=0.5)), col="black", lwd=3, lty=3)
lines(predict(loess(Carp3[,2]~Carp3[,5],span=0.5)), col="black", lwd=3, lty=4)
lines(predict(loess(Carp4[,2]~Carp4[,5],span=0.5)), col="black", lwd=3, lty=5)

legend(0,1,legend=c("Selfed w. Mort.", "Selfed no Mort.", "Full Sib", "Half Sib", "First Cousin", 
	"Minimum"), lwd=c(3,3,3,3,3,NA), lty=c(1,1,3,4,5,NA), pch=c(NA,NA,NA,NA,NA,5), 
	col=c("black","gray50","black","black","black","black"), bty="n", cex=1.2)

axis(1, at=c(0,20,40,60,80,100), labels=F)
axis(2, at=c(0,0.2,0.4,0.6,0.8,1.0), labels=c(0,0.2,0.4,0.6,0.8,1.0))

mtext("Carpinteria", side=2, line=5, cex=1.5)

#Goleta
par(mar=c(1,7,0,0), bty="n")

plot(predict(loess(Gol1[,2]~Gol1[,5],span=0.5)), col="black", type="l", lwd=3, lty=1, xlab="",
	ylab="Absolute Difference", ylim=c(0,1.0), xlim=c(0,100), axes=F)
lines(which.min(predict(loess(Gol1[,2]~Gol1[,5],span=0.5))), min(predict(loess(Gol1[,2]~Gol1[,5],span=0.5))),
	type="p", col="black", pch=5, cex=1.5)
lines(predict(loess(Gol0[,2]~Gol0[,5],span=0.5)), col="gray50", lwd=3, lty=1)
lines(predict(loess(Gol2[,2]~Gol2[,5],span=0.5)), col="black", lwd=3, lty=3)
lines(predict(loess(Gol3[,2]~Gol3[,5],span=0.5)), col="black", lwd=3, lty=4)
lines(predict(loess(Gol4[,2]~Gol4[,5],span=0.5)), col="black", lwd=3, lty=5)

axis(1, at=c(0,20,40,60,80,100), labels=F)
axis(2, at=c(0,0.2,0.4,0.6,0.8,1.0), labels=c(0,0.2,0.4,0.6,0.8,1.0))

mtext("Goleta", side=2, line=5, cex=1.5)


#Mohawk
par(mar=c(5,7,0,0), bty="n")

plot(predict(loess(Moh1[,2]~Moh1[,5],span=0.5)), col="black", type="l", lwd=3, lty=1, xlab="Percentage Non-Outcrossed", 
	ylab="Absolute Difference", ylim=c(0,1.0), xlim=c(0,100), axes=F)
lines(which.min(predict(loess(Moh1[,2]~Moh1[,5],span=0.5))), min(predict(loess(Moh1[,2]~Moh1[,5],span=0.5))),
	type="p", col="black", pch=5, cex=1.5)
lines(predict(loess(Moh0[,2]~Moh0[,5],span=0.5)), col="gray50", lwd=3, lty=1)
lines(predict(loess(Moh2[,2]~Moh2[,5],span=0.5)), col="black", lwd=3, lty=3)
lines(predict(loess(Moh3[,2]~Moh3[,5],span=0.5)), col="black", lwd=3, lty=4)
lines(predict(loess(Moh4[,2]~Moh4[,5],span=0.5)), col="black", lwd=3, lty=5)
	
axis(1, at=c(0,20,40,60,80,100), labels=c(0,20,40,60,80,100))
axis(2, at=c(0,0.2,0.4,0.6,0.8,1.0), labels=c(0,0.2,0.4,0.6,0.8,1.0))

mtext("Mohawk", side=2, line=5, cex=1.5)








#############################   PLOTS WITH CHI SQUARED GRAPHS   #####################################################

par(mfrow=c(3,2), bty="n")

### CARPINTERIA ###
#Selfing Plot, Smoothed Fit
plot(predict(loess(Carp1[,2]~Carp1[,5],span=0.5)), col="black", type="l", lwd=2, lty=1, xlab="",
	ylab="Absolute Difference in Proportion", main="a. Carpinteria", ylim=c(0,1.0), xlim=c(0,100))
#Lowest Value
lines(which.min(predict(loess(Carp1[,2]~Carp1[,5],span=0.5))), min(predict(loess(Carp1[,2]~Carp1[,5],span=0.5))),
	type="p", col="blue", pch=15)

#Selfing No Mortality Plots
lines(predict(loess(Carp0[,2]~Carp0[,5],span=0.5)), col="black", lwd=1, lty=1)

#Full Sib Plots	
	lines(predict(loess(Carp2[,2]~Carp2[,5],span=0.5)), col="blue", lty=6)
	
#Half Sib Plots
	lines(predict(loess(Carp3[,2]~Carp3[,5],span=0.5)), col="red", lty=2)
	
#First Cousin Plots	
	lines(predict(loess(Carp4[,2]~Carp4[,5],span=0.5)), col="green", lty=4)
	
# Legend WITH Selfing No Mortality Line
legend(0,1,legend=c("Selfed w. Mort.", "Selfed no Mort.", "Full Sib", "Half Sib", "First Cousin", 
	"Minimum"), lwd=c(2,1,1,1,1,NA), lty=c(1,1,6,2,4,NA), pch=c(NA,NA,NA,NA,NA,15), 
	col=c("black","black","blue","red","green","blue"), bty="n")


#ChiSq Plots
plot(0:100, rep(0.05,101), type="l", lty=3, lwd=2, col="red", ylim=c(0,1), xlab="",
	ylab="Chi Square Probability")
lines(predict(loess(Carp1[,4]~Carp1[,5],span=0.4)), lty=1, lwd=2)
lines(predict(loess(Carp2[,4]~Carp2[,5],span=0.4)), lty=6)

legend(60,1.0, legend=c("Selfed w. Mort.", "Selfed No Mort.", "All Other", "0.05 Critical Value"), lwd=c(2,1,1,2),
	lty=c(1,1,6,3), col=c("black","black","black","red"), bty="n")


### GOLETA ###
#Selfing Plot, Smoothed Fit
plot(predict(loess(Gol1[,2]~Gol1[,5],span=0.5)), col="black", type="l", lwd=2, lty=1, xlab="",
	ylab="Absolute Difference in Proportion", main="b. Goleta", ylim=c(0,1.0), xlim=c(0,100))
#Lowest Value
lines(which.min(predict(loess(Gol1[,2]~Gol1[,5],span=0.5))), min(predict(loess(Gol1[,2]~Gol1[,5],span=0.5))),
	type="p", col="blue", pch=15)

#Selfing No Mortality Plots
lines(predict(loess(Gol0[,2]~Gol0[,5],span=0.5)), col="black", lwd=1, lty=1)

#Full Sib Plots	
	lines(predict(loess(Gol2[,2]~Gol2[,5],span=0.5)), col="blue", lty=6)
	
#Half Sib Plots
	lines(predict(loess(Gol3[,2]~Gol3[,5],span=0.5)), col="red", lty=2)
	
#First Cousin Plots	
	lines(predict(loess(Gol4[,2]~Gol4[,5],span=0.5)), col="green", lty=4)
	

#ChiSq Plots
plot(0:100, rep(0.05,101), type="l", lty=3, lwd=2, col="red", ylim=c(0,1), xlab="",
	ylab="Chi Square Probability")
lines(predict(loess(Gol1[,4]~Gol1[,5],span=0.3)), lty=1, lwd=2)
lines(predict(loess(Gol2[,4]~Gol2[,5],span=0.3)), lty=6)


### MOHAWK ###
#Selfing Plot, Smoothed Fit
plot(predict(loess(Moh1[,2]~Moh1[,5],span=0.5)), col="black", type="l", lwd=2, lty=1, xlab="Percentage Non-Outcrossed", 
	ylab="Absolute Difference in Proportion", main="c. Mohawk", ylim=c(0,1.0), xlim=c(0,100))
#Lowest Value
lines(which.min(predict(loess(Moh1[,2]~Moh1[,5],span=0.5))), min(predict(loess(Moh1[,2]~Moh1[,5],span=0.5))),
	type="p", col="blue", pch=15)

#Selfing No Mortality Plots
lines(predict(loess(Moh0[,2]~Moh0[,5],span=0.5)), col="black", lwd=1, lty=1)

#Full Sib Plots	
	lines(predict(loess(Moh2[,2]~Moh2[,5],span=0.5)), col="blue", lty=6)
	
#Half Sib Plots
	lines(predict(loess(Moh3[,2]~Moh3[,5],span=0.5)), col="red", lty=2)
	
#First Cousin Plots	
	lines(predict(loess(Moh4[,2]~Moh4[,5],span=0.5)), col="green", lty=4)
	

#ChiSq Plots
plot(0:100, rep(0.05,101), type="l", lty=3, lwd=2, col="red", ylim=c(0,1), xlab="Percentage Non-Outcrossed", 
	ylab="Chi Square Probability")
lines(predict(loess(Moh1[,4]~Moh1[,5],span=0.3)), lty=1, lwd=2)
lines(predict(loess(Moh0[,4]~Moh0[,5],span=0.3)), lty=1, lwd=1)
lines(predict(loess(Moh2[,4]~Moh2[,5],span=0.3)), lty=6)

#which(Moh0[,4]>=0.05)