#Compute the expected change in frequency of a risky versus certain reproductive strategy
#assuming either fixed population size or variable population size

#generate number of offspring produced by parents with "risky" genotype
risky <- function(parents = 8){2*rbinom(1, parents, 0.5)}


deltap <- 

fixmut <- function(s = 0.01, N = 10,...) {
	#compute and store all post-selection frequencies
	p.star <- rep(0,2*N)
	for(j in 1:(2*N-1)){
		p.star[j] <- j*(1+s)/(j*s + 2*N)
		}
	
	#generate entries of the matrix Q
	#use the R function that computes the binomial probability density
	Q <- matrix(0, nrow = (2*N-1), ncol = (2*N-1))
	for(i in 1:(2*N-1)){
		for(j in 1:(2*N-1)) {
			Q[i, j] <- dbinom(j, 2*N, p.star[i])
			}
		}
	
	#generate entries fo the vector r (probs. of fixation in one step)
	r <- rep(0,(2*N-1))
	for(i in 1:(2*N-1)) {
		r[i] <- dbinom(2*N, 2*N, p.star[i])
		}
	
	#compute u, vector of ultimate fixation probabilities
	u <- rep(0,(2*N-1))
	I <- diag(2*N-1) #the identity matrix
	u <- solve(I - Q, r)
	u[1] #output the fixation probability of a new, unique mutation
	}

################################

############## 2N = 4 ##########
#neutral mutation
neut <- fixmut(s = 0, N = 2)  #should be 1/(2*N)
neut

#advantageous mutation with s = 0.01
fixmut(s = 0.01, N = 2)
fixmut(s = 0.01, N = 2)/neut #fixation prob relative to that of a neutral mutation

#deleterious mutation with s = -0.01
fixmut(s = -0.01, N = 2)
fixmut(s = -0.01, N = 2)/neut #fixation prob relative to that of a neutral mutation


############## 2N = 40 ##########
#neutral mutation
neut <- fixmut(s = 0, N = 20)  #should be 1/(2*N)
neut

#advantageous mutation with s = 0.01
fixmut(s = 0.01, N = 20)
fixmut(s = 0.01, N = 20)/neut #fixation prob relative to that of a neutral mutation

#deleterious mutation with s = -0.01
fixmut(s = -0.01, N = 20)
fixmut(s = -0.01, N = 20)/neut #fixation prob relative to that of a neutral mutation


############## 2N = 400 ##########
#neutral mutation
neut <- fixmut(s = 0, N = 200)  #should be 1/(2*N)
neut

#advantageous mutation with s = 0.01
fixmut(s = 0.01, N = 200)
fixmut(s = 0.01, N = 200)/neut #fixation prob relative to that of a neutral mutation

#deleterious mutation with s = -0.01
fixmut(s = -0.01, N = 200)
fixmut(s = -0.01, N = 200)/neut #fixation prob relative to that of a neutral mutation