# input parameters
p=0.1     	# selected proportion
varam=0.1 	# genetic variancve in phenotype
varav=0.05 	# genetic variance Ve
varp=1.0   	# phenotypic variance
rgamav=0.0   	# genetic correlation between am and av
ar=0.25		# additive genetic relationship between selection candidate and the group of relatives (0.25 for halfsib and 0.5 for parent-offspring)
aw=0.25		# additive genetic relationship among the relatives in the family
n=50		# number of relatives

# derived parameters
Ve=varp-varam
h2v=varav/(2*(varp^2)+3*varav) #heritability Ve
tv=aw*h2v	# intraclass correlation for variance

# accuracy of selection with relatives (question 6)
rih<-ar*sqrt(h2v)*sqrt(n/(1+(n-1)*tv))

# in a loop
n=0
accuracies=data.frame()
for (i in 1:20) {
n=n+10
rih<-ar*sqrt(h2v)*sqrt(n/(1+(n-1)*tv))
accuracies <- rbind(accuracies,data.frame(Nhalfsibs=n,accuracy=rih ))
}
directory=getwd()
write("OUTPUT: accuracies written to file","accuracies.txt",append=FALSE)
write(c(directory,date()),"accuracies.txt",append=TRUE)
write.table(accuracies,"accuracies.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File accuracies.txt written to directory ",directory))


# response in variance (question 7, 8)
n=50
rih<-ar*sqrt(h2v)*sqrt(n/(1+(n-1)*tv))
x<-qnorm(p,mean=0,sd=1,lower.tail=FALSE,log.p=FALSE) # truncation point normal distribution
d<-dnorm(x,mean=0,sd=1,log=FALSE) #value of the normal density function at truncation point x
i<-d/p # selection intensity
delta_Av<-i*rih*sqrt(varav)
Ve_1<-Ve+delta_Av



# accuracy, response in Ve and heritability Ve in the next 5 generations after selection in generation 0
results <- data.frame()
for (gen in 1:5) {
 varp1<-varam+Ve_1
 h2v=varav/(2*(varp1^2)+3*varav) #heritability Ve
 tv=aw*h2v	# intraclass correlation for variance
 rih<-ar*sqrt(h2v)*sqrt(n/(1+(n-1)*tv))
 delta_Av<-i*rih*sqrt(varav)

 results <- rbind(results,data.frame(generation=gen,accuracy=rih,h2v=h2v,Ve_1,varp1 ))
 Ve_1<-Ve_1+delta_Av
}

directory=getwd()
write("OUTPUT: variances after after selection","results.txt",append=FALSE)
write(c(directory,date()),"results.txt",append=TRUE)
write.table(results,"results.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File results.txt written to directory ",directory))

########################multi-trait selection###########

# multi-trait selection: response to selection in mean and variance with sibs/progeny (question 4)
# input pai*rameters
p=0.1     	# selected proportion
varam=0.1 	# genetic variancve in phenotype
varav=0.05 	# genetic variance Ve
varp=1.0   	# phenotypic variance
rgamav=0.0   	# genetic correlation between am and av
ar=0.5		# additive genetic relationship between selection candidate and the group of relatives (0.25 for halfsib and 0.5 for parent-offspring)
aw=0.25		# additive genetic relationship among the relatives in the family
n=100	     	# number of relatives
v_am<-1.0	# economic value Am
v_av<--1.0	# economic value Av

# derived parameters
Ve=varp-varam
h2v=varav/(2*(varp^2)+3*varav) #heritability Ve
tv=aw*h2v	# intraclass correlation for variance
h2m=varam/varp
tm=aw*h2m

x<-qnorm(p,mean=0,sd=1,lower.tail=FALSE,log.p=FALSE) # truncation point normal distribution
d<-dnorm(x,mean=0,sd=1,log=FALSE) #value of the normal density function at truncation point x
i<-d/p # selection intensity


P<-array(0,c(2,2))
P[1,1]<-((1+(n-1)*tm)*varp)/n
P[1,2]<-((3+aw*(n-3))*rgamav*sqrt(varav*varam))/n
P[2,1]<-P[1,2]
P[2,2]<-((1+(n-1)*tv)*(2*varp^2+3*varav))/n

invP<-solve(P)

G<-array(0,c(2,2))
G[1,1]<-ar*varam
G[1,2]<-ar*rgamav*sqrt(varav*varam)
G[2,1]<-ar*rgamav*sqrt(varav*varam)
G[2,2]<-ar*varav

v<-array(0,c(2))
v[1]<-v_am
v[2]<-v_av

b<-array(0,c(2))

b<-invP%*%G%*%v

response<-array(0,c(2))
sigmaI<-sqrt(t(b)%*%P%*%b)

response<-i*(t(b)%*%G)/sigmaI[1,1]
# genetic response is i*b*G/standard deviation index; first value is response in phenotype; second value is response in Ve
response

# question 10 when changing rgamav
rgamav=-1
results_rg=data.frame()
for (ii in 1:19) {
rgamav=rgamav+0.1
P<-array(0,c(2,2))
P[1,1]<-((1+(n-1)*tm)*varp)/n
P[1,2]<-((3+aw*(n-3))*rgamav*sqrt(varav*varam))/n
P[2,1]<-P[1,2]
P[2,2]<-((1+(n-1)*tv)*(2*varp^2+3*varav))/n

invP<-solve(P)

G<-array(0,c(2,2))
G[1,1]<-ar*varam
G[1,2]<-ar*rgamav*sqrt(varav*varam)
G[2,1]<-ar*rgamav*sqrt(varav*varam)
G[2,2]<-ar*varav

v<-array(0,c(2))
v[1]<-v_am
v[2]<-v_av

b<-array(0,c(2))

b<-invP%*%G%*%v

response<-array(0,c(2))
sigmaI<-sqrt(t(b)%*%P%*%b)

response<-i*(t(b)%*%G)/sigmaI[1,1]

results_rg <- rbind(results_rg,data.frame(rgamav=rgamav,response_am=response[1],response_av=response[2] ))
}
directory=getwd()
write("OUTPUT: results_rg written to file","results_rg.txt",append=FALSE)
write(c(directory,date()),"results_rg.txt",append=TRUE)
write.table(results_rg,"results_rg.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File results_rg.txt written to directory ",directory))

# question 11 when changing rgamav
v_av=1
rgamav=0.5
results_ev=data.frame()
for (ii in 1:101) {
v_av=v_av-1
P<-array(0,c(2,2))
P[1,1]<-((1+(n-1)*tm)*varp)/n
P[1,2]<-((3+aw*(n-3))*rgamav*sqrt(varav*varam))/n
P[2,1]<-P[1,2]
P[2,2]<-((1+(n-1)*tv)*(2*varp^2+3*varav))/n

invP<-solve(P)

G<-array(0,c(2,2))
G[1,1]<-ar*varam
G[1,2]<-ar*rgamav*sqrt(varav*varam)
G[2,1]<-ar*rgamav*sqrt(varav*varam)
G[2,2]<-ar*varav

v<-array(0,c(2))
v[1]<-v_am
v[2]<-v_av

b<-array(0,c(2))

b<-invP%*%G%*%v

response<-array(0,c(2))
sigmaI<-sqrt(t(b)%*%P%*%b)

response<-i*(t(b)%*%G)/sigmaI[1,1]

results_ev <- rbind(results_ev,data.frame(v_av=v_av,response_am=response[1],response_av=response[2] ))
}
directory=getwd()
write("OUTPUT: results_ev written to file","results_ev.txt",append=FALSE)
write(c(directory,date()),"results_ev.txt",append=TRUE)
write.table(results_ev,"results_ev.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File results_ev.txt written to directory ",directory))